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Range  Expansion  or  Range  Shift?  Population  Genetics  and  Historic 
Range  Data  Analyses  of  the  Predatory  Benthic  Sea  Slug 
Phidiana  hiltoni  (Mollusca,  Gastropoda,  Nudihranchia) 
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^Department  of  Biological  Sciences,  California  State  Polytechnic  University,  3801  West 

Temple  Avenue,  Pomona,  CA  91768 

^Department  of  Ecology  and  Evolutionary  Biology,  University  of  California,  Los 

Angeles,  CA  90095 

^Marine  Science  Institute,  University  of  California,  Santa  Barbara,  CA  93106 
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Abstract. — Phidiana  hiltoni  is  a  conspicuous  nudibranch  sea  slug  native  to  the  north¬ 
eastern  Pacific  Ocean.  Over  the  past  thirty  years  the  range  of  P.  hiltoni  has  expanded 
about  200  km  northward,  but  the  mechanism  that  facilitated  this  expansion  is  poorly 
understood.  In  this  study,  we  use  mtDNA  and  microsatellite  data  to  investigate  the 
population  structure  of  P.  hiltoni  in  its  historical  range  as  well  as  in  recently  colonized 
localities.  Microsatellite  analyses  reveal  little  to  no  genetic  structure  and  thus  high  gene 
flow  throughout  the  range  of  P.  hiltoni.  This  is  consistent  with  mtDNA  analysis  results, 
which  revealed  shared  haplotypes  between  Southern,  Central  and  Northern  popula¬ 
tions.  However,  AMOVA  of  mtDNA  data  did  recover  some  genetic  structure  among 
geographic  regions.  This,  along  with  same  group  memberships  in  the  microsatellite 
data  of  individuals  from  sites  like  Cave  Landing,  suggest  a  certain  degree  of  local 
recruitment  and  reduced  vagility.  Recently  established  populations  in  Northern  Cal¬ 
ifornia  contain  two  unique  mtDNA  haplotypes  that  are  not  present  elsewhere,  but 
microsatellite  data  do  not  differentiate  these  from  other  populations.  The  mismatch 
between  mtDNA  and  microsatellite  data  could  be  explained  by  the  mating  system  of 
this  aggressive,  hermaphroditic  species  as  well  as  the  sporadic  nature  of  the  northward 
dispersal.  Analyses  of  historical  abundance  data  of  P.  hiltoni  suggest  a  population  de¬ 
cline  in  Southern  California.  Together,  these  results  suggest  a  northward  population 
shift,  rather  than  a  range  expansion,  possibly  related  to  ongoing  changes  in  nearshore 
oceanographic  conditions  in  the  region. 


Rising  ocean  temperatures  driven  by  global  climate  change  are  having  dramatic  impacts 
on  coastal  ecosystems  around  the  world  (McGowan  et  al.  1998;  Sorte  et  al.  2011).  One 
of  the  most  noticeable  effects  is  the  poleward  range  expansion  of  certain  species  (Dawson 
et  al.  2010;  Sorte  et  al.  201 1;  Sunday  et  al.  2012;  Canning-Clode  and  Carlton  2017).  Par¬ 
ticularly  problematic  are  range  expansions  of  predatory  species,  which  can  have  significant 
impacts  on  the  trophic  structure  of  newly  colonized  ecosystems  (Zeidberg  and  Robinson 
2007;  Gallardo  et  al.  2016).  However,  not  all  of  these  range  expansions  are  permanent; 
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some  result  from  regular  oscillations  in  ocean  temperatures  (e.g.,  El  Nino  events).  In  these 
cases,  populations  often  return  to  their  original  range  following  ephemeral  warming  events, 
making  it  difficult  to  attribute  in  the  short  term  any  particular  range  shift  to  longer  term 
climate  change  (Schultz  et  al.  201 1).  Poleward  range  expansions  can  also  be  accompanied 
by  extirpation  at  lower  latitudes,  resulting  in  shifts  at  both  ends  of  species  ranges  (Parme¬ 
san  et  al.  1999;  Bates  et  al.  2014).  Range  shifts  may  constitute  a  more  pervasive  indication 
of  permanent  changes  in  the  ecological  structure  of  biotas  as  they  can  be  more  difficult  to 
reverse  (Parmesan  et  al.  1999;  Schultz  et  al.  2011).  However,  range  shifts  are  difficult  to 
detect  and  precisely  quantify,  particularly  in  marine  species  with  low  abundance  and/or 
population  densities  (Bates  et  al.  2015). 

Phidiana  hiltoni  is  a  relatively  large  and  conspicuous  aeolid  nudibranch  native  to 
the  northeastern  Pacific  Ocean.  Like  most  nudibranchs,  P.  hiltoni  is  a  simultaneous 
hermaphrodite,  but  the  mating  behavior  of  this  species  is  poorly  understood.  The  diet 
of  P.  hiltoni  consists  mostly  of  hydroids  and  other  cnidarians;  however,  individuals  of 
this  species  are  known  to  attack  and  consume  other  sea  slugs,  particularly  small,  soft- 
bodied  aeolids  and  dendronotaceans,  including  conspecifics  (Goddard  et  al.  2011).  Phidi- 
anci  hiltoni  has  relatively  large  eggs  and  lecithotrophic  larval  development  (Goddard  2004); 
its  larvae  do  not  need  to  feed  in  the  water  column  and  are  capable  of  settlement  and  meta¬ 
morphosis  within  a  day  or  two  of  hatching.  Thus,  compared  to  planktotrophic  species, 
dispersal  by  the  larvae  of  P.  hiltoni  is  greatly  reduced.  Historically,  P.  hiltoni  was  found  as 
far  south  as  Isla  Cedros  off  the  coast  of  Baja  California,  Mexico  and  as  far  north  as  Pacific 
Grove,  California  (Goddard  et  al.  2011).  In  1977,  Phidiana  hiltoni  was  discovered  north  of 
Monterey  Bay  (Goddard  et  al.  201 1).  Once  across  Monterey  Bay,  P  hiltoni  rapidly  made 
its  way  up  the  coast;  it  was  found  just  north  of  San  Francisco  Bay  (Duxbury  Reef,  Marin 
County)  in  1992  and  now  is  present  as  far  north  as  Bodega  Bay,  California,  representing  a 
200-km  northward  range  expansion  in  40  years  (Goddard  et  al.  2011;  Goddard  et  al.  2018). 
At  Duxbury  Reef,  P.  hiltoni  quickly  became  the  dominant  sea  slug,  with  apparent  negative 
impacts  on  other  nudibranch  species,  likely  through  a  combination  of  direct  predation  and 
competition  for  shared  hydroid  prey  (Goddard  et  al.  2011). 

The  mechanism  behind  the  range  expansion  of  P.  hiltoni  is  not  well  understood,  but 
has  been  potentially  linked  to  warming  coastal  waters  and  shifts  in  ocean  currents  along 
the  California  coast  (Schultz  et  al.  201 1).  Changes  in  ocean  circulation,  which  drives  lar¬ 
val  transport,  can  potentially  increase  the  risk  of  species  introductions  and/or  dispersals 
(Harley  et  al.  2006;  Sorte  et  al.  2011;  Wilson  et  al.  2016).  Whereas  the  dispersal  potential  of 
species  with  planktonic  feeding  larvae  is  relatively  well  understood  (Scheltema  1986),  less  is 
known  about  how  species  with  lecithotrophic  development  may  respond  to  oceanographic 
changes.  A  majority  of  lecithotrophic  sea  slug  species  are  found  in  warmer,  nutrient-poor 
waters,  where  having  a  short-lived,  non-feeding  larval  stage  can  lead  to  reductions  in  larval 
mortality  at  the  cost  of  reduced  fecundity  and  vagility  (Goddard  2004;  Goddard  and  Her- 
mosillo  2008).  Phidiana  hiltoni  is  one  of  the  few  lecithotrophically  developing  nudibranchs 
found  in  temperate  waters  in  the  Northeast  Pacific  Ocean  (Goddard  2004),  and  its  recent 
range  expansion  may  reflect  long-term  changes  in  regional  nearshore  circulation  regimes 
and  productivity  (Rebstock  2003).  These  factors  make  P.  hiltoni  a  particularly  interesting 
system  for  studying  the  complex  interactions  between  climate  change,  range  shifts,  and 
marine  invasion  biology. 

If  the  dispersal  of  P.  hiltoni  has  been  facilitated  by  changes  of  oceanographic  regimes 
and  warming  waters  at  the  northern  edge  of  its  range,  the  most  likely  source  for  the  new 
populations  north  of  Monterey  Bay  are  Central  California  populations.  However,  it  is  also 
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possible  that  individuals  from  father  south  were  introduced  into  northern  California  either 
by  larval  dispersal  or  human  activities.  The  lecithotrophic  larval  development  of  P.  hiltoni 
makes  it  an  ideal  candidate  for  ballast  water  dispersal,  and  less  likely  to  disperse  long 
distances  naturally  in  response  to  environmental  changes.  Two  of  the  busiest  commercial 
ports  in  North  America  are  located  in  California,  with  the  Los  Angeles-Long  Beach  Har¬ 
bor  well  within  the  historic  range  of  P.  hiltoni,  and  the  Oakland-San  Francisco  Harbor 
(San  Francisco  Bay)  in  the  center  of  the  extended  portion  of  the  range.  However,  the  ab¬ 
sence  of  P  hiltoni  from  San  Francisco  Bay  along  with  its  prevalence  in  open-coast  rocky 
reefs  (Goddard  et  al.  2011)  contradicts  the  ballast  water  introduction  hypothesis.  Another 
possible  vector  for  the  spread  of  P.  hiltoni  could  be  small  vessel  traffic  between  regional 
ports  and  harbors  (Wasson  et  al.  2001),  but  the  mobile  hunting  behavior  of  this  species 
makes  it  an  unlikely  fouling  organism.  Regardless  of  the  mechanism  of  dispersal,  Schultz 
et  al.  (201 1 )  noted  that  P  hiltoni  has  persisted  at  higher  latitudes  despite  ocean  temperature 
fluctuations  from  El  Niiio/La  Nina  cycles,  suggesting  this  species  may  be  an  indicator  of 
faunal  range  shifts  due  to  climate  change. 

Although  the  range  expansion  and  some  of  the  associated  ecological  effects  of  P.  hiltoni 
are  well  documented  (Goddard  et  al.  201 1),  many  questions  remain.  No  genetic  studies 
have  been  conducted  on  P.  hiltoni,  thus  the  population  structure  of  the  species  is  unknown, 
hampering  our  ability  to  understand  the  mechanisms  of  dispersal.  Also,  very  little  atten¬ 
tion  has  been  paid  to  the  southern  range  limit  of  P.  hiltoni,  leaving  unanswered  the  question 
as  to  whether  recent  observations  indicate  a  northern  range  expansion  or  overall  range 
shift.  In  the  present  study,  we  examine  population  structure  in  P.  hiltoni,  explore  the  ge¬ 
netic  signature  of  its  range  expansion,  and  hypothesize  possible  dispersal  mechanisms  into 
Northern  California.  Additionally,  we  reviewed  historical  collection  data  near  the  south¬ 
ern  range  of  P.  hiltoni  in  order  to  document  population  density  changes  that  may  help  to 
understand  the  population  dynamics  of  this  species. 


Materials  and  Methods 

Whole  specimens  and  tissue  samples  of  Phicliana  hiltoni  (Table  SI)  were  obtained  from 
different  sources  and  various  locations  along  the  California  coast  (Fig.  1).  Some  individu¬ 
als  were  collected  at  the  shoreline  during  low  tide,  other  specimens  and  tissue  samples  were 
provided  by  colleagues  or  obtained  from  collections  of  the  Natural  History  Museum  of  Los 
Angeles  County  (LACM),  the  Santa  Barbara  Museum  of  Natural  History  (SBMNH)  and 
the  California  Academy  of  Sciences  (CASIZ).  Fieldwork  was  conducted  under  the  Cal¬ 
ifornia  Department  of  Fish  and  Wildlife  permit  #13256.  Specimens  collected  in  the  field 
were  preserved  in  95%  ethanol  and  deposited  at  the  California  State  Polytechnic  University 
Invertebrate  Collection  (CPIC). 

DNA  was  extracted  from  sixty  specimens  (Table  1)  using  a  DNeasy  Blood  and  Tis¬ 
sue  Kit  (Qiagen,  Valencia,  CA)  using  standard  protocols  provided  by  the  manufacturer. 
A  fragment  of  the  cytochrome  c  oxidase  subunit  I  (COI)  mitochondrial  gene  was  am¬ 
plified  and  sequenced  using  universal  primers  (LCO1490  5'-GGTCAACAAATCATAAA 
GATATTGG-3',  HC02198  5'-TAAACTTCAGGGTGACCAAAAAATCA-3')  (Folmer 
et  al.  1994).  PCR  reaction  conditions  were  as  follows:  denaturation  at  95°C  for  3  min, 
followed  by  35  cycles  of  94°C  for  45  sec,  50°C  for  45  sec,  and  72°C  for  2  min  with  fi¬ 
nal  elongation  at  72°C  for  10  min.  Successful  DNA  amplification  was  confirmed  using  an 
agarose  gel  electrophoresis  and  ethidium  bromide.  The  PCR  products  were  purified  with  a 
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Fig.  1.  Map  of  the  historic  and  expanded  range  of  Phidiana  liiltoiii  along  the  western  coast  of  North 
America.  Sampling  localities  are  indicated  by  arrows  and  coded  with  different  grey  tones  by  geographic 
regions. 


GeneJET  PCR  Purification  Kit  (Fermentas,  Waltham,  MA)  using  standard  protocols  and 
were  sent  for  sequencing  to  Source  Bioscience  Inc.  (Santa  Fe  Springs,  CA). 

Sequences  were  assembled  and  aligned  using  Geneious  v8.1.8  (Kearse  et  al.  2012).  The 
geographic  distribution  of  mtDNA  haplotypes  was  visualized  by  producing  a  haplotype 
network  using  the  program  PopArt  vl.7  (Feigh  and  Bryant  2015)  using  the  TCS  option. 
Haplotypes  were  pattern-coded  by  locality.  Genetic  structure  within  and  among  popula¬ 
tions  and  among  groups  was  examined  using  analysis  of  molecular  variance  (AMOVA)  as 
implemented  in  Arlequin  v3.5  (Excoffier  and  Fischer  2010).  Three  different  AMOVA  anal¬ 
yses  were  run  to  test  for  the  effects  of  arranging  populations  into  different  groups  based 
on  the  distribution  of  genetic  variation.  In  the  first  AMOVA  test,  populations  established 
after  1977  were  included  in  the  Northern  California  group  while  historic  range  populations 
were  divided  into  two  groups:  Central  California  (populations  north  of  Point  Conception) 
and  Southern  California  (populations  south  of  Point  Conception);  this  is  the  hypothesized 
biogeographic  structure  if  Point  Conception  acts  as  a  biogeographic  barrier  (Blanchette 
et  al.  2008).  To  examine  possible  genetic  similarities  between  newly  formed  populations 
in  Northern  California  and  those  from  southern  Monterey  Bay,  two  additional  AMOVAs 
were  run  with  different  group  arrangements,  to  examine  whether  this  resulted  in  different 
distributions  of  genetic  diversity  among  groups  and  among  populations  within  groups.  In 
the  second  AMOVA,  the  groups  were  kept  the  same  except  central  California  populations 
were  split  into  two  groups.  Northern  Central  California  (southern  Monterey  Bay)  and 
Southern  Central  California  (populations  further  south),  see  Blanchette  et  al.  (2008).  In 
the  third  AMOVA  populations  from  southern  Monterey  Bay  were  pooled  with  Northern 
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California.  Significance  of  the  AMOVAs  was  tested  using  16,000  permutations  of  individ¬ 
uals  between  groups.  Arlequin  v3.5  was  also  used  to  calculate  pairwise  Ost  between  pop¬ 
ulations  (1,000  permutations).  Because  the  Southern  California  populations  (Naples  and 
San  Clemente)  were  represented  by  one  sequence  each,  populations  were  pooled  together. 

Microsatellite  loci  were  identified  through  sequencing  by  synthesis  with  a  MiSeq  plat¬ 
form  (Illumina,  Inc.,  San  Diego,  CA).  DNA  was  extracted  from  a  single  specimen  collected 
from  Pillar  Point,  California  (CASIZ  190249),  tagged  with  a  unique  barcode  during  library 
preparation,  and  pooled  with  other  samples  for  Illumina  sequencing.  Sequencing  was  con¬ 
ducted  at  the  UCLA  Genotyping  and  Sequencing  Core  facility.  Automated  screening  of 
sequences  for  tetranucleotide  repeats  and  primer  design  were  performed  simultaneously  in 
MSATCOMMANDER  vl.0.8  (Faircloth  2008).  Twenty-two  primer  pairs  were  purchased 
from  Eurofins  (Louisville,  KY)  with  a  M13  tail  added  to  the  3'  end  of  each  forward  primer 
sequence. 

Five  of  the  twenty-two  primer  pairs  were  tested  with  ten  specimens  that  consistently  am¬ 
plified  for  mtDNA  to  determine  the  PCR  protocol.  PCR  protocol  settings  for  the  primers 
were  optimized  from  a  standard  protocol  by  adjusting  the  annealing  temperatures  and 
elongation  times  until  amplification  was  achieved.  The  PCR  Master  Mix  for  each  locus  in 
these  tests  included  the  forward  primer  with  a  M13  tail,  reverse  primer,  and  BSA  (bovine 
serum  albumin)  and  used  Thermo  Fisher  Platinum  Hot  Start  PCR  Master  Mix.  The  final 
PCR  reaction  conditions  were  as  follows:  denaturation  at  95°C  for  15  min,  followed  by 
30  cycles  of  94°C  for  30  sec,  60-65°C  for  30  sec,  and  68°C  for  30  sec  with  final  elonga¬ 
tion  at  60°C  for  10  min.  All  twenty-two  primer  pairs  were  tested  with  ten  specimens  that 
consistently  amplified  for  mtDNA  using  the  above  conditions.  Of  the  twenty-two  primer 
pairs  tested,  ten  polymorphic  loci  amplified  reliably.  Using  the  ten  reliable  primer  pairs 
(Table  S2)  and  the  above  amplification  conditions,  PCR  was  carried  out  with  fifty-seven 
specimens.  The  PCR  Master  Mix  for  each  locus  in  this  final  round  now  included  a  fluores¬ 
cent  M 13  tag  (5'-[6-FAM]  AGGGTTTTCCCAGTCACGACGTT-3')  along  with  the  orig¬ 
inal  components.  Genotyping  was  outsourced  to  Laragen  Incorporated  (Culver  City,  CA). 
Genotypes  were  scored  using  the  Microsatellite  Analysis  External  Plugin  v  1.4.4  imple¬ 
mented  in  Geneious  v8.1.8  using  the  Two  Surrounding  Peaks  setting  (Kearse  et  al.  2012). 

In  total,  fifty-five  individual  specimens  were  genotyped  for  all  10  microsatellite  loci.  This 
is  a  small  sample  size  for  this  type  of  study,  but  specimens  were  difficult  to  obtain  in  the 
field.  Collecting  sea  slugs  is  serendipitous  in  nature  and  after  two  years  of  fieldwork  only  a 
small  number  of  specimens  was  obtained.  Additionally,  most  museum  specimens  examined 
were  unsuitable  for  molecular  work.  Population  subdivision  in  the  nuclear  genome  was  in¬ 
ferred  using  STRUCTURE  v2.3.4  (Pritchard  et  al.  2000)  with  the  default  parameters;  5 
replicates  for  each  value  of  K  were  run  for  1,000,000  MCMC  iterations  following  a  burn-in 
period  of  100,000.  To  detect  the  true  number  of  clusters  {K)  using  the  Evanno  Method 
(Evanno  et  al.  2005)  the  result  file  from  STRUCTURE  was  processed  with  STRUCTURE 
Harvester  vO.6.9.84  (Earl  and  vonHoldt  2012).  Using  the  selected  K  value  (3)  the  result¬ 
ing  files  were  processed  with  CLUMPP  vl.1.2  (Rosenberg  et  al.  2002)  and  Distruct  vl.l 
(Rosenberg,  2004)  to  generate  a  graphic  display  of  the  population  structure.  AMOVA  and 
Fst  pairwise  genetic  differentiation  comparisons  between  populations  were  conducted  fol¬ 
lowing  the  same  methodology  as  in  the  mtDNA  analyses.  Microsatellite  data  were  also 
analyzed  via  Discriminant  Analysis  of  Principal  Components  (DAPC)  using  the  adegent 
package  in  R  (Jombart  et  al.  2010). 

To  determine  whether  the  abundance  of  P.  hiltoni  in  Southern  California  has  changed 
since  the  mid-20^*^  century,  counts  of  nudibranchs  by  James  R.  Lance  dating  from  1953 
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to  2001'  at  six  rocky  intertidal  sites  in  San  Diego  County  (Point  Loma,  Hill  Street, 
False  Point,  Bird  Rock,  Windansea,  and  South  Casa  Reef)  were  examined  and  ana¬ 
lyzed.  The  data  for  P.  hiltoni  were  extracted  and  grouped  by  site,  decade,  and  before 
and  after  1963,  the  year  when,  excepting  one  brief  trip  to  Bird  Rock  in  1956,  Lance 
started  sampling  outer  coast  sites  in  San  Diego  County  other  than  Point  Loma.  Counts 
made  on  consecutive  or  near-consecutive  dates  at  any  given  site  were  excluded  from 
analysis  in  order  to  reduce  autocorrelation  in  the  data;  the  count  retained  was  the  one 
with  the  highest  number  of  P.  hiltoni.  Twenty-one  additional  counts  conducted  by  ei¬ 
ther  JG  or  CK  from  2000  to  2016  at  4  of  the  same  sites  (Point  Loma,  Hill  Street, 
Bird  Rock,  and  South  Casa  ReeQ  were  also  included  in  the  analysis.  A  Wilcoxon 
sign-rank  test  was  implemented  in  JMP  vl3,  SAS  Institute  Inc.  (Cary,  NC)  and  used 
to  compare  the  number  of  P.  hiltoni  found  at  Point  Loma  before  and  after  1963. 
Additional  information  on  the  recent  occurrence  of  P.  hiltoni  in  San  Diego  County 
was  obtained  from  the  website  iNaturalist  (https://www.inaturalist.org/taxa/48724- 
Phidiana-hiltoni)  and  the  species  database  on  Divebums,  a  San  Diego  dive  website 
(http:// species.divebums.com/index. php?l=sciname&n=Phidiana%20hiltoni),  and  con¬ 
firmed  by  the  authors. 


Results 

The  haplotype  network  of  the  mitochondrial  COI  gene  recovered  five  distinct  haplo- 
types  (Fig.  2).  Thirty-three  individuals  spanning  all  seven  populations  share  the  most  com¬ 
mon  haplotype.  Three  haplotypes  diverge  from  the  most  common  haplotype  by  only  two 
nucleotides.  The  most  common  of  these  three  haplotypes  was  found  exclusively  in  four¬ 
teen  individuals  from  Pillar  Point,  Northern  California.  The  other  two  haplotypes  are  only 
found  in  specimens  originating  from  Cayucos  and  Cave  Landing,  Central  California.  An 
additional  specimen  from  Cave  Landing  possessed  a  haplotype  diverging  from  the  most 
common  haplotype  by  three  nucleotides. 

In  the  first  AMOVA  test  (Northern  California:  Pillar  Point;  Central  California:  Cave 
Landing,  Jalama  Beach,  Cayucos,  and  Carmel  Point)  most  of  the  genetic  variation  is  re¬ 
covered  among  groups  (70.05%)  and  within  populations  (30.59%),  with  virtually  no  vari¬ 
ation  among  populations  within  groups  (-0.63%)  (Table  2).  In  the  second  AMOVA  test 
(Northern  California:  Pillar  Point;  Northern  Central  California:  Carmel  Point;  Southern 
Central  California:  Jalama  Beach,  Cave  Landing,  Cayucos)  most  of  the  genetic  variation 
is  again  among  groups  (65.94%)  and  within  populations  (32.82%)  and  very  little  variation 
among  populations  within  groups  (1.24%)  (Table  2).  In  the  third  AMOVA  test  (Northern 
California:  Pillar  Point,  Carmel  Point;  Central  California:  Jalama  Beach,  Cave  Landing, 
Cayucos)  most  of  the  genetic  variation  is  found  again  among  groups  (50.47%)  and  within 
populations  (33.5%)  however  genetic  variation  is  found  among  populations  within  groups 
(16.04%)  (Table  2). 

A  pairwise  Ost  test  was  run  on  all  populations  and  resulted  in  relatively  high  values  be¬ 
tween  Pillar  Point  and  all  other  populations  and  also  between  Carmel  Point  and  Cayucos. 
The  only  significant  difference  in  genetic  variation  found  was  between  Pillar  Point,  North¬ 
ern  California  and  each  of  the  four  Central  California  populations:  Jalama  Beach,  Cave 


'  Goddard,  JH.R.  2013.  Opisthobranch  gastropods  observed  on  the  outer  coast  of  San  Diego  County, 
California  by  James  R.  Lance,  1953-2001.  knb.298.2.  [online]  California  Academy  of  Sciences.  Available 
https://knb.ecoinformatics.Org/knb/metacat/knb.298.2/knb  [2017  Jim  15], 
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Fig.  2.  Haplotype  network  of  COI  mitochondrial  sequences  generated  with  PopArt  vl.7.  Each  circle 
represents  a  unique  haplotype  and  its  area  is  proportional  to  the  number  of  specimens  sequenced  with  that 
haplotype.  Each  pattern  represents  the  geographic  origin  of  the  individual  specimens,  as  indicated  in  the 
legend.  Isolate  codes  are  indicated  next  to  each  haplotype. 


Landing,  Cayucos  and  Carmel  Point  (Ost  =  0.69, =  0.03;  <5>st  =  0.75,/?  =  0.00;  cE>st  = 
0.69,  p  =  0.00;  OsT  =  0.67,  p  =  0.04  respectively)  (Table  3).  This  suggests  genetic  differen¬ 
tiation  between  Pillar  Point  and  the  other  populations.  However,  this  result,  as  well  as  the 
lack  of  significant  differentiation  among  other  pairwise  comparisons,  must  be  interpreted 
with  caution  due  to  the  limited  number  of  samples  from  all  collection  sites  except  for  Pillar 
Point  and  Cave  Landing. 

Analysis  of  microsatellite  data  with  Structure  Harvester  using  Evanno’s  method 
(Evanno  et  al.  2005),  a  maximum  value  of  the  rate  of  change  (A AT)  in  the  log  probabil¬ 
ity  of  data  was  obtained  at  AT  =  3  (Fig.  3B).  These  three  recovered  clusters  are  unevenly 
distributed  among  geographic  regions  with  no  obvious  geographic  subdivision  (cluster  1 : 
red,  cluster  2:  blue,  cluster  3:  yellow,  Fig.  3A).  Moreover,  all  individuals  exhibit  a  non-zero 
probability  of  belonging  to  any  one  of  the  three  clusters. 

All  AMOVA  tests  with  different  groupings  produced  very  similar  results;  the  overwhelm¬ 
ing  majority  of  the  genetic  variation  was  recovered  within  populations  (92.32-93.61%)  and 
some  among  populations  within  groups  (6.51-9.98%),  with  virtually  no  variation  among 
groups  (-3.59-1.17%)  (Table  4).  Pairwise  Fst  comparisons  produced  very  low  values 
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Table  3.  Ost  pairwise  comparison  values  for  mitochondrial  haplotype  data  obtained  with  Arlequin  v3.5 
(lower  triangular)  and  associated  p  values  (upper  triangular).  Significant  values  (p  <  0.05)  in  bold. 


Pillar 

Point 

Carmel 

Point 

Cayucos 

Cave 

Landing 

Jalama 

Beach 

Naples 

San 

Clemente 

Pillar  Point 

— 

0.03062 

0.00019 

0.0000 

0.03537 

0.09994 

0.99994 

Carmel  Point 

0.69490 

— 

0.21358 

0.32000 

0.99994 

0.99994 

0.99994 

Cayucos 

0.75345 

0.62791 

— 

0.99994 

0.99994 

0.99994 

0.99994 

Cave  Landing 

0.68730 

0.29687 

-0.05381 

- 

0.99994 

0.99994 

0.99994 

Jalama  Beach 

0.66771 

0.0000 

0.0000 

-0.33043 

— 

0.99994 

0.99994 

Naples 

0.61778 

0.0000 

0.0000 

0.0000 

0.0000 

- 

0.99994 

San  Clemente 

0.61778 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

- 

across  the  entire  range  suggesting  little  to  no  genetic  dififerentiation  between  populations 
(Table  5).  DAPC,  which  attempts  to  group  individuals  using  a  k-means  clustering  algo¬ 
rithm,  suggests  that  the  entire  metapopulation  of  P.  hiltoni  cannot  be  divided  into  more 
than  one  group  based  on  microsatellite  data  (Fig.  4). 

From  1953  to  1962,  Jim  Lance  sampled  for  nudibranchs  on  the  outer  coast  of  San 
Diego  County  only  at  Point  Loma,  where  he  found  at  least  one  Phidicma  hiltoni  dur¬ 
ing  half  of  his  28  trips  during  that  period  (Fig.  5A).  In  1964  he  began  to  sample  addi¬ 
tional  outer  coast  sites.  Since  then,  only  two  more  P.  hiltoni  were  found  at  Point  Loma 
(Figure  5A  &  B),  one  by  Lance  in  July  1968  and  one  by  JG  in  June  2001.  Similarly, 
P.  hiltoni  was  observed  on  only  about  10%  of  the  trips  to  each  of  the  other  five  outer 
coast  sites  (Fig.  5 A),  and  was  found  in  lower  numbers  per  trip  than  had  been  seen  in 
the  earlier  period  at  Point  Loma  (Fig.  5B).  Significantly  fewer  P.  hiltoni  were  found  at 
Point  Loma  after  1963  than  before  (Wilcoxon  rank  sum  test,  P  =  0.002)  (Figs.  5B,  C). 
Except  for  the  1990s,  when  only  10  total  trips  were  made  to  two  sites,  P.  hiltoni  has 
been  found  in  San  Diego  County  in  low  numbers  in  each  of  the  decades  since  the  1960s 
(Fig.  5C).  Finally,  P.  hiltoni  has  been  photographed  subtidally  in  San  Diego  County 
at  least  10  times  since  2005  (http://species.divebums.com/index.php?!  =  sciname&n  = 
Phidiana%20hiltoni;  https://www.inaturalist.org/taxa/48724-Phidiana-hiltoni). 


Fig.  3.  Genetic  clustering  analysis  for  the  entire  data  set  of  57  individuals  and  10  microsatellite  re¬ 
gions  as  estimated  by  STRUCTURE  v2.3.4.  A.  Genetic  clustering  plot  for  K  =  3  clusters,  generated  with 
CLUMPP  vl.1.2.  Each  grey  tone  represents  a  difierent  genetic  cluster.  Bar  graphs  show  average  posterior 
probability  of  membership  (^v-axis)  of  each  individual.  Populations  are  delimited  by  dark  vertical  lines.  B. 
Graph  of  AK  =  mean  ( |  V'iK)  \ )  /  sdlLIA'))  as  a  function  of  K  (potential  number  of  genetic  clusters)  gen¬ 
erated  by  STRUCTURE  Harvester  vO.6.9.84.  The  most  likely  number  of  clusters  is  indicated  by  the  modal 
value,  in  this  case  K=  3. 


Table  4.  AMOVA  test  results  for  microsatellite  genotype  data  with  three  separate  population  groupings  obtained  with  Arlequin  v3.5,  significant  values  {p  <  0.05)  in 
bold.  Grouping  1:  Northern  California  (Pillar  Point),  Central  California  (Jalama  Beach,  Cave  Landing,  Cayucos,  Carmel  Point),  and  Southern  California  (San  Clemente, 
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Table  5.  Fst  pairwise  comparison  values  for  microsatellite  genotype  data,  obtained  with  Arlequin  v3.5 
(lower  triangular)  and  associated  p  values  (upper  triangular),  significant  values  (/;  <  0.05)  in  bold.  The  two 
southern  California  populations  of  Naples  and  San  Clemente  were  combined  into  one  due  to  low  sampling 
numbers. 


Pillar  Point 

Scott  Creek 

Carmel  Point 

Cave  Landing 

Jalama  Beach 

S.  California 

Pillar  Point 

— 

0.93776 

0.01968 

0.64388 

0.39049 

0.35181 

Scott  Creek 

0.13279 

— 

0.99994 

0.99994 

0.99994 

0.99994 

Carmel  Point 

0.14426 

0.35183 

- 

0.00306 

0.00444 

0.10229 

Cave  Landing 

0.01779 

0.17272 

0.15117 

- 

0.06255 

0.05986 

Jalama  Beach 

0.04219 

0.23112 

0.19336 

0.05260 

— 

0.19821 

S.  California 

0.12656 

0.36283 

0.30272 

0.15087 

0.15965 

- 

Discussion 

The  California  coast  is  a  prime  example  of  a  region  where  climate  change  is  impact¬ 
ing  native  marine  ecosystems  through  changes  in  ocean  temperatures,  seawater  chemistry, 
and  coastal  current  regimes  (Barry  et  al.  1995;  Sagarin  et  al.  1999;  Harley  et  al.  2006). 
With  ocean  temperatures  increasing,  species  are  predicted  to  shift  their  ranges  poleward,  a 
trend  that  has  already  been  observed  across  a  wide  range  of  taxonomic  groups,  mostly  at 
temperate  latitudes  (Dawson  et  al.  2010;  VanDerWal  et  al.  2012).  It  is  difficult  to  predict 
how  environmental  change,  including  biotic  exchanges  resulting  from  species  range  shifts 
and  introductions,  will  affect  ecological  systems.  Heilman  et  al.  (2008)  emphasized  that 


Fig.  4.  Results  of  the  DAPC  analysis  indicating  that  the  data,  when  analyzed  as  principal  components, 
cannot  be  divided  into  more  than  one  group. 
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A  0,6 


Point  Loma  Hill  Street  False  Point  Bird  Rock  Windansea  So.  Casa  Reef 


C  3,5  -1 


3.0  - 
2.5  - 


1950’s 


1960’s 


■  Point  Loma  »  Other  5  sites 


1970’s  1980’s  1990’s  2000’s  201 0’s 


Fig.  5.  Occurrence  of  Phidiana  hiltoni  at  six  rocky  intertidal  sites  in  San  Diego  County,  1953-2016.  A. 
Frequency  of  occurrence  by  site,  1953-1962  and  1964-2016.  Number  of  trips  to  each  site  shown  above  bars. 
B.  Number  of  P.  hiltoni  found  on  each  trip,  1953-2016.  C.  Mean  number  (±  SE)  oi  P.  hiltoni  found  per  trip 
at  Point  Loma  and  other  sites  combined,  by  decade. 
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global  climate  change  compounds  this  difficulty  because  it  transforms  transport  and  in¬ 
troduction  mechanisms,  impacts  the  distribution  of  existing  invasive  species,  and  alters  the 
effectiveness  of  control  strategies.  Efforts  to  manage  and  conserve  marine  ecosystems  in  the 
face  of  climate  change  will  require  improvements  to  the  existing  predictive  framework  to 
aid  in  preventing  future  introductions  (Harley  et  al.  2006).  In  this  context,  understanding 
the  mechanisms  behind  range  expansions  (or  shifts)  of  individual  species  will  contribute 
to  a  larger  body  of  evidence,  critically  important  for  predicting  the  biological  effects  of 
climate  change. 

The  data  presented  in  this  paper  provide  insight  into  the  processes  underlying  the  range 
expansion  in  P.  hiltoui.  The  analysis  of  mitochondrial  DNA  (mtDNA)  sequence  data  of 
Phidicma  hiltoni  across  both  its  historical  range  and  extended  range  in  Northern  California 
revealed  limited  genetic  structure.  The  COI  haplotype  network  shows  very  little  polymor¬ 
phism  with  only  five  haplotypes  in  total  being  recovered  (Fig.  2).  The  most  common  of 
these  haplotypes  is  shared  among  individuals  from  all  seven  populations,  including  the 
newly  formed  populations  in  Northern  California.  However,  there  are  a  few  haplotypes 
only  detected  in  certain  populations.  For  example,  one  of  the  two  haplotypes  found  at  Pil¬ 
lar  Point,  north  of  Monterey  Bay,  was  found  nowhere  else.  Because  populations  north  of 
Monterey  Bay  did  not  exist  prior  to  1977  (Goddard  et  al.  2011),  it  is  likely  that  this  unique 
haplotype  is  also  present  south  of  Monterey  Bay  but  at  such  low  frequency  that  it  has  not 
been  yet  detected,  and  may  have  become  more  common  in  the  recently  colonized  popu¬ 
lations  due  to  founder  effects.  To  better  understand  the  geographic  structure  of  P.  hiltoni 
based  on  mtDNA,  three  AMOVA  analyses  were  run,  each  with  different  population  group¬ 
ings.  In  each  of  the  three  groupings,  the  highest  percent  variation  was  consistently  found 
among  groups  (Table  2).  However,  there  is  a  decrease  in  this  percent  variation  as  the  group¬ 
ings  structure  are  altered  (70.05%,  65.94%,  50.47%  respectively),  suggesting  that  the  first 
grouping  arrangement  (in  which  the  Northern  California  group  includes  only  populations 
found  in  the  extended  range)  best  represents  population  genetic  structure  according  to 
mtDNA  data.  Pairwise  Ost  comparisons  agree  with  the  geographic  structure  recovered  in 
the  AMOVAs  (Table  3).  Relatively  high  Ost  values  indicating  genetic  differentiation  were 
found  between  Pillar  Point  and  the  Central  and  Southern  California  populations,  and  be¬ 
tween  Carmel  Point  and  the  two  populations  of  Cayucos  and  Cave  Fanding.  This  is  con¬ 
sistent  with  the  haplotype  network  results,  showing  that  Pillar  Point,  Cayucos  and  Cave 
Fanding  possess  divergent  haplotypes  from  the  most  common  haplotype  found  across  the 
range  of  P.  hiltoni. 

While  analyses  of  mtDNA  sequences  suggest  population  structure  and  genetic  differen¬ 
tiation  among  groups,  no  genetic  pattern  corresponding  to  geography  was  detected  using 
microsatellite  data.  Across  the  range  of  P.  hiltoni.  Structure  analyses  showed  several  indi¬ 
viduals  have  nearly  identical  probabilities  of  cluster  membership,  a  pattern  that  is  partic¬ 
ularly  apparent  at  Cave  Fanding.  This  genetic  uniformity  could  be  an  indication  of  self¬ 
recruitment  at  Cave  Fanding,  where  the  concavity  in  the  coastline  (accentuated  by  the 
750-meter  long  rock  jetty  on  Point  San  Fuis)  may  encourage  larval  retention.  Goddard 
et  al.  (2011)  suggested  that  population  structure  in  P.  hiltoni  should  be  affected  by  up- 
welling  shadows,  resulting  in  local  retention  of  short-lived  larvae,  particularly  at  the  north¬ 
ern  end  of  bights  along  the  coastline  (Graham  and  Fargier  1997;  Roughan  et  al.  2005). 
AMOVAs  on  microsatellite  data  included  the  same  three  distinct  groupings  as  mitochon¬ 
drial  analyses  except  for  the  addition  of  Scott  Creek  to  the  Northern  California  group  and 
the  removal  of  Cayucos  (where  microsatellite  data  were  not  successfully  recovered).  In  all 
AMOVA  tests,  most  genetic  variance  was  found  within  populations  (93.61%,  92.32%  and 
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93.  1 1%,  respectively)  (Table  4).  Consistently  low  variation  among  groups  is  indicative  of 
high  levels  of  gene  flow  across  the  range,  as  seen  in  the  mtDNA  data.  Pairwise  compar¬ 
isons  of  Fst  values  from  microsatellite  data  are  also  consistent  with  high  gene  flow  across 
the  range,  with  generally  low  Fst  values  between  populations  (Table  5).  The  slightly  higher 
Fst  values  found  between  Scott  Creek  and  Carmel  Point,  as  well  as  between  Carmel  Point 
and  Southern  California  are  unreliable  due  to  the  low  sample  sizes  from  Scott  Creek  and 
Southern  California. 

AMOVA  and  pairwise  <J>st  comparisons  using  mtDNA  data  suggest  Northern  Califor¬ 
nia  populations  are  genetically  distinct,  and  consequently  the  origin  of  the  recently  founded 
populations  remains  unclear.  On  the  contrary,  microsatellite  data  indicate  very  high  levels 
of  gene  flow  in  P.  hiltoni,  with  Northern  California  populations  genetically  indistinguish¬ 
able  from  those  in  Central  California.  This  discrepancy  of  results  from  mtDNA  vs.  nuclear 
data  might  be  explained  by  the  reproductive  behavior  of  P.  hiltoni.  Rutowksi  (1983)  found 
that  species  with  high  rates  of  cannibalism  require  several  couplings,  usually  with  differ¬ 
ent  mates,  in  order  for  all  the  eggs  in  the  egg  mass  to  be  fertilized.  It  is  very  likely  that 
P.  hiltoni  (considering  the  pugnacious  and  cannibalistic  nature  of  this  species)  exhibits  ab¬ 
breviated  coupling  times,  requiring  several  mates  to  fully  fertilize  egg  masses.  Additionally, 
P.  hiltoni  exhibits  locally  high  population  densities  in  Central  California,  potentially  facili¬ 
tating  promiscuity  and  thus  intense  sperm  competition.  The  main  consequence  of  this  mat¬ 
ing  system  would  be  that  egg  masses  produced  by  a  single  individual  (functional  female) 
will  be  sired  by  multiple  partners,  generating  offspring  with  genetically  identical  mitochon¬ 
dria  but  different  nuclear  alleles.  If  recently  established  populations  in  Northern  California 
are  the  result  of  sporadic  events  involving  dispersal  of  larvae  (or  rafting  of  egg  masses) 
produced  by  a  small  number  of  females,  it  is  likely  that  founder  effects  facilitated  retention 
of  these  rare  mitochondrial  haplotypes.  While  the  diversity  of  nuclear  alleles  should  also 
decrease  due  to  drift  and  founder  effects,  the  polyandrous  mating  system  of  P.  hiltoni  has 
the  potential  to  mitigate  these  effects.  Under  this  scenario,  Northern  California  popula¬ 
tions  could  harbor  mitochondrial  haplotypes  that  were  previously  very  rare  in  the  historic 
range,  while  nuclear  alleles  from  the  source  population  are  more  broadly  represented.  If 
this  hypothesis  is  correct,  additional  sampling  from  Central  California  should  detect  all 
or  most  Northern  California  haplotypes.  An  alternative  explanation  is  that  the  lack  of  ge¬ 
netic  structure  in  microsatellite  data  is  an  artifact  of  the  limited  sample  size.  If  this  is  the 
case,  additional  sampling  across  the  range  of  P.  hiltoni  would  improve  the  reliability  of  the 
results  of  this  study. 

Another  outstanding  question  is  what  mechanism(s)  allowed  P.  hiltoni  to  cross  Monterey 
Bay  starting  in  the  late  1970s,  or  what  physical  or  biological  barriers  restricted  the  prior 
range  of  this  species.  If  recently  established  populations  in  Northern  California  are  indeed 
the  result  of  sporadic  dispersal  by  a  limited  number  of  individuals,  this  would  suggest  that 
this  dispersal  was  not  triggered  by  a  gradual  process,  such  as  increasing  ocean  temper¬ 
atures,  but  instead  by  the  temporary  or  intermittent  opening  of  a  corridor.  One  distinct 
possibility  is  that  a  weakening  of  the  upwelling  shadow  in  Monterey  Bay  (Graham  and 
Largier  1997)  due  to  climate  change  may  have  facilitated  this  process.  Pennington  et  al. 
(2000)  documented  decadal-scale  changes  in  the  oceanographic  conditions  near  the  center 
of  Monterey  Bay  region  consistent  with  those  described  for  the  1976-77  climate  shift  in 
the  North  Pacific  Ocean.  These  include  increased  stratification  of  surface  waters,  warmer, 
less  productive  waters  during  non-upwelling  seasons  and  a  later  onset  of  upwelling.  These 
changes  may  have  reduced  larval  retention  in  northern  Monterey  Bay,  allowing  P.  hiltoni 
to  disperse  into  Northern  California. 
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Few  studies  have  examined  the  population  genetics  of  other  marine  invertebrate  taxa 
that  have  experienced  recent  range  expansions  or  range  shifts.  In  California,  a  similar  study 
by  Dawson  et  al.  (2010)  examined  three  hypotheses/scenarios  that  explain  the  causes  of 
range  limits  of  species,  and  concluded  that  in  the  volcano  barnacle,  Tetraclitci  ruhescens,  the 
northern  range  boundary  is  maintained  by  migration  load  arising  from  flow  of  maladapted 
alleles  into  peripheral  locations.  Dawson  et  al.  (2010)  proposed  that  in  this  species  (with 
planktonic-dispersing  larvae),  environmental  amelioration,  likely  due  to  climate  change, 
resulted  in  a  reduction  of  the  strength  of  selection  against  immigrant  phenotypes  in  the 
northern  range  boundary,  allowing  the  species  to  expand  northward.  The  case  of  R  hiltoui 
is  very  different  in  several  respects,  but  the  main  difference  is  that  whereas  T.  ruhescens  was 
declining  near  its  northern  range  limit,  P.  hiltoni  has  been  and  remains  common.  Framing 
of  our  data  in  the  three  scenarios/hypotheses  proposed  by  Dawson  et  al.  (2010)  suggest 
physical  barriers  to  dispersal  is  the  most  likely  mechanism  that  historically  restricted  mi¬ 
gration  in  P.  hiltoni. 

This  study  included  a  sample  of  57  individual  specimens  collected  across  the  range  of 
P  hiltoni,  but  the  sample  size  from  Southern  California,  south  of  Point  Conception,  is 
small.  Only  two  individuals  were  collected  despite  a  substantial  collecting  effort  by  the  se¬ 
nior  author  in  this  region  during  two  consecutive  years.  This  suggests  that  P.  hiltoni  could 
have  become  rare  in  the  southern  portion  of  its  range.  However,  this  assumption  should 
be  interpreted  with  caution.  Bates  et  al.  (2015)  shown  that  abundance-related  species  de¬ 
tectability,  particularly  important  in  uncommon,  difficult-to-detect  marine  species  such  as 
P.  hiltoni,  has  the  potential  to  confound  our  understanding  of  the  true  location  of  range 
edges.  Bates  et  al.  (2015)  emphasized  the  importance  of  simulation  and  modeling,  but  also 
long-term  monitoring  with  consistent  sampling  effort  through  time.  In  this  case,  we  ana¬ 
lyzed  high-quality,  long-term  observational  data  collected  at  fairly  regular  intervals  from 
the  same  region.  These  data  suggest  that  the  historical  abundance  of  P.  hiltoni  in  San  Diego 
County,  especially  Point  Loma  (the  type  locality  of  P.  pngncix  Lance  1961  [  =  P.  hiltoni]} 
has  declined  (Figs.  4A-C).  Reasons  for  this  decline  remain  unknown  and  warrant  fur¬ 
ther  investigation,  although  rising  ocean  temperatures  appears  to  be  a  viable  hypothesis. 
Notably,  the  abundance  of  another  species  of  sea  slug,  Feliniare  cciliforniensis,  was  once- 
common  in  Southern  California,  but  became  extinct  there  in  the  1980’s  (Goddard  et  al. 
2013).  Although  individuals  of  F.  cciliforniensis  reappeared  in  2003  and  the  species  has 
since  been  found  in  a  few  isolated  localities  in  Southern  California  (Goddard  et  al.  2013; 
Hoover  2015),  its  populations  have  not  completely  recovered.  It  is  unclear  whether  there  is 
a  link  between  the  decline  of  these  two  ecologically  distinct  species,  but  if  there  is,  it  may 
be  a  symptom  of  larger  and  more  pervasive  environmental  change.  The  apparent  decline 
of  P.  hiltoni  in  Southern  California  along  with  its  dispersal  northward  needs  to  be  sub¬ 
stantiated  with  further  monitoring  and  additional  data  analyses  (Bates  et  al.  2015),  but  if 
confirmed,  would  suggest  this  is  a  true  poleward  range  expansion  rather  than  a  temporary 
shift  (Parmesan  et  al.  1999).  Understanding  the  process  by  which  P  ////rw?/ migrated  north¬ 
ward  may  provide  insight  as  to  how  other  benthic  organisms  will  respond  to  rising  ocean 
temperatures  and  changes  in  ocean  current  systems  (McGowan  et  al.  1998). 
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Appendix 


Table  SI .  Complete  list  of  specimens  sequenced  for  this  study,  including  isolate  number,  locality,  collec¬ 
tion  date,  and  GenBank  accession  numbers. 


Isolate 

Locality 

Collection  date 

GenBank  accession  # 

CKll 

San  Clemente  Island,  CA 

1961 

MK333330 

CK29 

Cave  Landing,  CA 

10/8/2014 

MK333291 

CK30 

Cave  Landing,  CA 

10/8/2014 

MK333292 

CK31 

Cave  Landing,  CA 

10/8/2014 

MK333293 

CK32 

Cave  Landing,  CA 

10/8/2014 

MK333294 

CK33 

Cave  Landing,  CA 

10/8/2014 

MK333295 

CK34 

Naples,  CA 

12/2009 

MK333313 

CK46 

Tarantula  Reef,  Jalama  Beach,  CA 

12/15/2009 

MK333331 

CK47 

Tarantula  Reef,  Jalama  Beach,  CA 

12/15/2009 

MK333332 

CK48 

Tarantula  Reef,  Jalama  Beach,  CA 

12/15/2009 

MK333333 

CK49 

Tarantula  Reef,  Jalama  Beach,  CA 

12/15/2009 

MK333334 

CK50 

Tarantula  Reef,  Jalama  Beach,  CA 

12/15/2009 

MK333335 

CK51 

Tarantula  Reef,  Jalama  Beach,  CA 

12/15/2009 

MK333336 

CK52 

Tarantula  Reef,  Jalama  Beach,  CA 

12/15/2009 

MK333337 

CK53 

Tarantula  Reef,  Jalama  Beach,  CA 

12/15/2009 

MK333338 

CK60 

Pillar  Point,  CA 

6/18/2015 

MK333314 

CK61 

Pillar  Point,  CA 

6/18/2015 

MK333315 

CK62 

Pillar  Point,  CA 

6/18/2015 

MK333316 

CK63 

Pillar  Point,  CA 

6/18/2015 

MK333317 

CK64 

Pillar  Point,  CA 

6/22/2015 

MK333318 

CK65 

Pillar  Point,  CA 

6/22/2015 

MK333319 

CK66 

Pillar  Point,  CA 

6/22/2015 

MK333320 

CK67 

Pillar  Point,  CA 

6/22/2015 

MK333321 

CK68 

Pillar  Point,  CA 

6/22/2015 

MK333322 

CK69 

Pillar  Point,  CA 

6/22/2015 

MK333323 

CK70 

Pillar  Point,  CA 

6/22/2015 

MK333324 

CK71 

Pillar  Point,  CA 

6/22/2015 

MK333325 

CK73 

Pillar  Point,  CA 

6/22/2015 

MK333326 

CK74 

Pillar  Point,  CA 

6/22/2015 

MK333327 

CK75 

Pillar  Point,  CA 

6/22/2015 

MK333328 

CK76 

Pillar  Point,  CA 

6/22/2015 

MK333329 

CK77 

Cave  Landing,  CA 

1/16/2016 

MK333296 

CK78 

Cave  Landing,  CA 

1/16/2016 

MK333297 

CK79 

Cave  Landing,  CA 

1/16/2016 

MK333298 

CK80 

Cave  Landing,  CA 

1/16/2016 

MK333299 

CK81 

Cave  Landing,  CA 

1/16/2016 

MK333300 

CK82 

Cave  Landing,  CA 

1/16/2016 

MK333301 

CK83 

Cave  Landing,  CA 

1/16/2016 

MK333302 

CK84 

Cave  Landing,  CA 

1/16/2016 

MK333303 

CK86 

Cave  Landing,  CA 

1/16/2016 

MK333304 

CK87 

Cave  Landing,  CA 

1/16/2016 

MK333305 

CK88 

Cave  Landing,  CA 

1/16/2016 

MK333306 

CK89 

Cave  Landing,  CA 

2/5/2016 

MK333307 

CK93 

Cave  Landing,  CA 

2/5/2016 

MK333308 

CK94 

Cave  Landing,  CA 

2/5/2016 

MK333309 

CK98 

Cave  Landing,  CA 

2/5/2016 

MK333310 

CK117 

Carmel  Pt,  Monterey  Bay  CA 

— 

MK333289 

CK118 

Carmel  Pt,  Monterey  Bay  CA 

— 

MK333290 

CK130 

Cayucos,  CA 

— 

MK33331  1 

CK131 

Cayucos,  CA 

— 

MK333312 
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Table  S2.  List  of  primer  pairs  (and  their  sequences)  used  to  amplify  polymorphic  microsatellite  loci  in 
P.  hiltoui.  Bolded  portion  of  forward  primers  indicate  M13  tail. 


Primer 

Sequence 

Phil760625F 

Phil760625R 

Phil792112F 

Phil792112R 

Phil820905F 

Phil820905R 

Phil928092F 

Phil928092R 

Phil98151F 

Phil98151R 

Phill09255F 

Phill09255R 

Phil585958F 

Phil585958R 

Phil918696F 

Phil918696R 

Phill21774F 

Phill21774R 

Phil315595F 

Phil315595R 

AGGGTTTTCCCAGTCACGACGTTAACGTCGTCATGGAATTCACAG 

GTTTATTAATGGCGGCGATGTGAC 

AGGGTTTTCCCAGTCACGACGTTAACCAATCGACGACAAGCTAAC 

GTTTGTCTCCGTGTTAAGTGTTGC 

AGGGTTTTCCCAGTCACGACGTTACATTACTCCACTCGACTCAGG 

GTTTAGTCTCGGTCCATGAATCAGG 

AGGGTTTTCCCAGTCACGACGTTGATTCTATGCCACACACCTTGG 

GTTTAATGTATCTGCTTCATCCGTGC 

AGGGTTTTCCCAGTCACGACGTTAGAGGAATAGTCGCGGAACTAC 

GTTTCATCATTGCGTCAGATGTCC 

AGGGTTTTCCCAGTCACGACGTTCACACGTTCATACACTCACCTG 

GTTTAACACCGAGACAAGACATGC 

AGGGTTTTCCCAGTCACGACGTTACTCTCTCACACCTGTCAAGTC 

GTTTCACCTCAGTACAGTCTCGTG 

AGGGTTTTCCCAGTCACGACGTTACTCTCTCACACCTGTCAAGTC 

GTTTCACCTCAGTACAGTCTCGTG 

AGGGTTTTCCCAGTCACGACGTTGTCAAGTGAATAAGACGGCGAG 

GTTTCTGCCTGCTATACATCCATCC 

AGGGTTTTCCCAGTCACGACGTTGTAACACAGTGTCCGTATGTGG 

GTTTATCATTCTACGTGCATGCTGTC 
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Status  of  the  Endangered  Indian  Knob  Mountainbahn  Eriodictyon 
altissimum  (Namaceae)  in  Central  Coastal  Cabfornia 

Christopher  P.  Kofron,**  Connie  Rutherford,^  Lisa  E,  Andreano," 

Michael  J.  Walgren,“  and  Heather  Schneider"^ 

^  U.S.  Fish  and  Wildlife  Service,  2493  Portola  Road,  Suite  B,  Ventura,  CA  93003 
-California  Department  of  Parks  and  Recreation,  750  Hearst  Castle  Road, 

San  Simeon,  CA  93452 

^ Santa  Barbara  Botanic  Garden,  1212  Mission  Canyon  Road,  Santa  Barbara, 

CA  93105 

Abstract. — Indian  Knob  Mountainbahn  Eriodictyon  altissimum  (Namaceae)  is  a 
shrub  endemic  to  western  San  Luis  Obispo  County  in  central  coastal  California,  and 
little  has  been  published  regarding  it.  The  species  was  listed  as  endangered  under  the 
California  Endangered  Species  Act  in  1979  and  the  U.S.  Endangered  Species  Act  in 
1995.  At  Federal  listing  in  1995,  Indian  Knob  mountainbalm  was  known  from  six 
occurrences,  two  of  which  were  in  protected  areas,  with  a  total  population  estimate 
of  <600  individuals.  As  of  2019,  Indian  Knob  mountainbalm  is  known  from  seven 
occurrences,  six  of  which  are  in  protected  areas  and  one  (the  largest)  mostly  in  a  pro¬ 
tected  area,  with  a  total  population  count  of  6,489+  individuals  in  2016.  Two  occur¬ 
rences  are  likely  extirpated.  Indian  Knob  mountainbalm  is  considered  a  fire-adapted 
chaparral  plant.  Reproduction  is  reported  to  be  primarily  vegetative  by  underground 
rhizomes,  and  it  is  specialized  for  substrates  with  physical  disturbances,  including: 
steep  rocky  slopes,  cliff  faces,  fallen  rock  debris,  sand  dunes  (shifting  sand),  road¬ 
sides,  old  graded  substrates  such  as  dirt/rock  roads,  the  talus  of  graded  substrates, 
and  trails.  We  report  the  species  grows  up  to  5.5  m  tall  and  at  98  to  263  m  eleva¬ 
tion.  In  consideration  of  the  life  history  traits  used  by  Anacker  et  al.  (2013)  for  rare 
plants  in  California,  Indian  Knob  mountainbalm  would  be  considered  highly  vul¬ 
nerable  to  climate  change.  Using  the  international  standards  of  lUCN,  Indian  Knob 
mountainbalm  meets  the  criteria  for  classification  as  endangered  including  the  follow¬ 
ing:  geographic  range,  fragmented;  extent  of  occurrence,  34  km-  (<100  km-);  area 
of  occupancy,  <2.3  km-  (<10  km^);  and  quality  of  habitat,  continuing  to  decline 
(dense  vegetation,  lack  of  recent  fire).  Coordinated  conservation  and  research  are 
needed  to  further  understand  the  species,  and  to  restore  and  maintain  the  five  extant 
occurrences. 


Indian  Knob  mountainbalm  Eriodictyon  altissimum  (Namaceae,  Luebert  et  al.  2016; 
Fig.  1)  is  a  shrub  endemic  to  western  San  Luis  Obispo  County  in  central  coastal  Cal¬ 
ifornia  (Fig.  2).  The  species  was  listed  as  endangered  under  the  California  Endangered 
Species  Act  in  1979  (CDFW  2016)  and  the  U.S.  Endangered  Species  Act  in  1995  (USFWS 
1994).  It  is  recognized  also  as  a  IB.l  rare  plant  (seriously  threatened)  by  the  California 


*  Corresponding  author:  chris_kofron@fws.gov 


21 


SOUTHERN  CALIFORNIA  ACADEMY  OF  SCIENCES 


O') 


Fig.  1.  Left:  A  tall  individual  of  Indian  Knob  mountainbalin  Eriodictyou  altissimimi  with  flowers  in 
Morro  Dunes  Ecological  Reserve  East  in  Los  Osos.  San  Luis  Obispo  County.  California  (occurrence  6).  20 
April  2016.  Right:  An  extraordinarily  tall  individual  of  Indian  Knob  mountainbalin  rising  above  the  ma¬ 
ture  chaparral  comprised  predominantly  of  Santa  Margarita  manzanita  Arclostaphylos  pilosula  at  Baron 
Canyon  Ranch  Estates  (previously  occurrence  7,  now  part  of  occurrence  5),  San  Luis  Obispo  County,  Cali¬ 
fornia.  27  July  2016.  The  foreground  is  recently  cleared  and  graded.  Brandon  Sanderson  of  California  De¬ 
partment  of  Fish  and  Wildlife  is  1.88  m  tall. By  extrapolation,  the  individual  of  Indian  Knob  mountainbalm 
is  determined  to  be  ~5.5  m  tall. 


Native  Plant  Society^  At  Federal  listing  in  1995,  Indian  Knob  mountainbalm  was  known 
from  six  occurrences,  two  of  which  were  in  protected  areas,  with  a  total  population  es¬ 
timate  of  <600  individuals  (USFWS  1994).  Five  of  the  occurrences  were  small  and  re¬ 
ported  to  each  comprise  <50  individuals,  and  the  largest  occurrence  on  the  private  Guidetti 
Ranch  was  reported  to  comprise  <350  individuals  (USFWS  1994).  Three  of  these  occur¬ 
rences  were  on  private  properties,  two  were  on  land  owned  by  the  State  of  California 


'  California  Native  Plant  Society.  2018.  Eriodiclyou  ciltissimimi.  In:  Inventory  of  rare  and  endangered 
plants  of  California.  Available  (July  2018):  http://www.rareplants.cnps.org/detail/623.html. 
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Fig.  2.  The  known  geographic  distribution  of  Indian  Knob  mountainbalm  Eriodictyou  cdtissinunn  in 
western  San  Luis  Obispo  County  in  central  coastal  California.  Each  occurrence  is  indicated  by  a  black  dot 
(extant)  or  black  circle  (likely  extirpated)  with  the  number  assigned  by  the  California  Natural  Diversity 
Database.  Occurrence  8  is  new. 


(Montana  de  Oro  State  Park),  and  one  was  on  land  owned  by  the  County  of  San  Luis 
Obispo  and  used  for  wastewater  treatment.  The  primary  threat  at  the  time  of  listing  was 
habitat  loss  as  a  result  of  potential  development  (USFWS  1994),  specifically  potential  sur¬ 
face  mining  at  the  largest  occurrence  on  Indian  Knob  mountain,  and  potential  residential 
development  at  two  small  occurrences  in  Los  Osos.  Our  purpose  is  to  review  and  enhance 
the  knowledge  of  Indian  Knob  mountainbalm,  in  particular  its  distribution,  ecology,  abun¬ 
dance,  threats,  management  and  conservation  status  in  2019.  Herein  we  present  the  first 
comprehensive  survey  and  census  of  the  species  throughout  its  geographic  range. 

Materials  and  Methods 

Little  has  been  published  regarding  Indian  Knob  mountainbalm  (but  see  Wells  1962; 
Carlquist  et  al.  1983;  Hannan  2012).  We  use  our  newly-collected  data,  along  with  the  sparse 
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published  literature,  information  in  the  California  Natural  Diversity  Database  (CDFG 
2006;  CDFW  2017),  unpublished  reports  and  other  documents  in  files  of  the  USFWS 
to  provide  a  status  review  of  the  species.  We  searched  for  Indian  Knob  mountainbalm, 
and  censused  and  mapped  the  known  occurrences  in  2016,  2017  and  2018.  Our  searches 
included  nearby  and  potential  habitat,  and  we  state  the  specific  dates  of  our  searches  in 
Appendix  1.  In  accordance  with  lUCN  (2014)  recommended  guidelines  for  colonial  and 
modular  organisms,  we  counted  the  numbers  of  stems  (including  ramets)  of  Indian  Knob 
mountainbalm  as  we  walked  along  the  trails  and  through  the  open  areas  at  each  known 
occurrence.  In  addition,  at  previous  occurrence  7  (Baron  Canyon  Ranch  Estates,  now  part 
of  occurrence  5),  we  used  binoculars  to  count  the  number  of  stems  in  the  distance  rising 
above  the  vegetation  canopy,  and  we  also  counted  the  number  of  stems  along  the  primary 
ridge  road  (Balm  Ridge  Way)  as  we  drove  slowly  on  it  (due  to  time  limitation).  In  order 
to  determine  potential  pollinators,  we  opportunistically  collected  insects  inside  the  corolla 
tubes  of  flowers  of  Indian  Knob  mountainbalm  during  April  and  May  2016.  We  measured 
the  stem  height  and  basal  diameter  of  several  extraordinarily  large  individuals  that  were 
accessible  to  us.  Using  all  available  information,  we  summarize  the  body  of  knowledge  of 
Indian  Knob  mountainbalm,  including  its  distribution,  known  occurrences,  ecology,  abun¬ 
dance,  threats,  management  and  conservation  status  in  2019.  We  considered  a  location  with 
the  species  as  a  separate  occurrence  if  it  was  >0.4  km  from  the  nearest  occurrence  (CDFG 
2011).  Latitude,  longitude  and  elevation  were  determined  with  a  global  positioning  sys¬ 
tem  (GPS)  device  in  the  field,  or  they  were  determined  using  Google  Earth  aerial  imagery. 
Common  and  Latin  names  of  plants  follow  Baldwin  et  al.  (2012).  Stated  areas  (ha)  of  prop¬ 
erties  are  from  records  of  the  County  of  San  Luis  Obispo.  In  Appendix  1 ,  we  discuss  each 
occurrence  in  sequence  according  to  its  assigned  number  in  the  California  Natural  Diver¬ 
sity  Database  (CDFG  2006;  CDFW  2017).  Thus,  there  may  be  no  numerical  sequence  for 
neighboring  occurrences.  Our  stated  numbers  of  plants  are  the  numbers  that  we  recorded 
but  they  should  be  considered  best  approximates  due  to  the  counting  methods  and  difficult 
landscapes,  and  in  this  context  we  refer  to  plants,  individuals  and  stems  interchangeably. 
The  findings  and  conclusions  in  this  article  are  those  of  the  authors  and  do  not  necessarily 
represent  the  views  of  the  USFWS. 


Results  and  Discussion 

Review  of  Indian  Knob  Mount  ainhahn. — Indian  Knob  mountainbalm  is  one  of  1 1  species 
(plus  additional  subspecies  and  varieties)  in  the  genus  Eviodietyon  that  occurs  in  southwest 
United  States  and  Mexico  (Hannan  2012).  It  is  a  diffusely-branched,  evergreen  shrub  with 
a  main  stem,  long  narrow  leaves  (5  to  9  cm  by  2  to  4  mm),  and  lavender  bell-shaped  flowers 
(11  to  16  mm  long;  USFWS  1998;  Hannan  2012;  Fig.  3).  The  species  possesses  under¬ 
ground  rhizomes  and  aboveground  stems.  Although  reported  heights  of  stems  are  up  to 
4  m  and  basal  stem  diameters  up  to  12  cm  (Wells  1962;  CDFG  2005;  Hannan  2012),  we 
observed  one  individual  that  was  ~5.5  m  tall  (Fig.  1)  and  another  individual  with  basal 
stem  diameter  of  ~20  cm.  Individuals  just  several  centimeters  tall  are  capable  of  produc¬ 
ing  flowers  (Fig.  3).  We  collected  a  total  of  56  insects  inside  the  corolla  tubes  of  flowers 
of  Indian  Knob  mountainbalm,  with  ants  (43%)  and  beetles  (23%)  comprising  66%  of  the 
sample.  The  insect  species  were  diverse  and  included  24  ants  (individuals),  13  beetles,  8 
bumblebees,  6  flies  (5  bee  flies  +  1  other  fly),  4  butterflies,  and  1  aphid.  These  data  indicate 
that  Indian  Knob  mountainbalm  receives  many  floral  visitors  but  it  is  unknown  which  are 
effective  pollinators. 
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Fig.  3.  Upper:  Flowers  of  Indian  Knob  mountainbalm  Eriodictyou  altissimimi  on  a  tall  individual, 
southwest  Guidetti  Ranch,  San  Luis  Obispo  County.  California  (occurrence  5).  4  May  2016.  Lower:  An 
individual  of  Indian  Knob  mountainbalm  that  is  only  several  centimeters  tall  and  with  flowers  at  southwest 
Guidetti  Ranch  (occurrence  5),  San  Luis  Obispo  County,  California,  5  May  2016. 
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Wells  (1962)  described  Indian  Knob  mountainbalm  as  rapidly  growing  and  short  lived, 
although  USFWS  (1998)  considered  the  species  long  lived  in  light  of  the  slow  growing 
lichens  on  some  individuals.  We  and  Dr.  Jon  Keeley  (U.S.  Geological  Survey)  examined  the 
growth  rings  in  the  dead  stems  of  several  large  individuals.  Despite  wood  rot  destruction, 
one  individual  had  >50  visible  growth  rings.  Although  we  do  not  know  if  one  growth  ring 
equates  to  1  yr,  J.  Keeley  (pers.  comm.  2018)  speculated  the  stems  were  >50  yr  of  age. 
Carlquist  and  Eckhart  (1984)  reported  that  the  underground  rhizomes  of  Eriodiclyon  are 
the  more  permanent  parts  of  the  plant,  and  they  may  be  older  than  the  aboveground  stems. 

Indian  Knob  mountainbalm  is  a  pioneer,  early  successional  or  edge  species  in  chaparral 
communities.  At  Indian  Knob,  Wells  (1962)  referred  to  it  as  a  pioneer  species  that  is  “ag¬ 
gressive  on  roadsides  with  numerous  young  plants  invading  such  disturbed  sites.”  CDFW 
(2017)  summarized  the  habitat  as  ridges  in  open  disturbed  areas  within  chaparral  on  sand¬ 
stone  and  shale,  and  openings  in  chaparral  on  stabilized  sand  dunes.  Specifically,  we  ob¬ 
served  the  species  on  sandstone  and  sand,  including  steep  rocky  slopes,  cliff  faces,  fallen 
rock  debris,  sand  dunes,  roadsides,  old  graded  substrates  such  as  dirt/rock  roads,  the  talus 
of  graded  substrates,  and  trails.  In  particular,  we  did  not  observe  the  species  on  shale  at 
any  occurrence. 

Although  Indian  Knob  mountainbalm  is  considered  to  be  fire-adapted  (Wells  1962; 
USFWS  1994,  1998,  2013b;  CDFG  2005),  there  is  no  direct  evidence  because  of  the 
absence  of  recent  fire  at  the  known  occurrences.  However,  in  general,  establishment  of 
seedlings  of  chaparral  plants  is  abundant  after  fire  but  uncommon  under  mature  chaparral 
(Keeley  1984).  Reproduction  of  Indian  Knob  mountainbalm  is  believed  to  be  primarily 
vegetative  by  underground  rhizomes  (USFWS  1994,  1998,  2013b;  CDFG  2005;  Fig.  4),  al¬ 
though  there  is  no  direct  evidence  to  support  this.  Because  little  is  published  about  Indian 
Knob  mountainbalm,  information  is  frequently  extrapolated  from  what  is  known  about 
other  species  in  the  genus,  such  as  Lompoc  yerba  santa  E.  capitatum  and  California  yerba 
Santa  E.  californicwn.  Lompoc  yerba  santa,  which  occurs  in  Santa  Barbara  County,  is 
self-incompatible,  meaning  that  pollen  from  genetically  different  plants  is  needed  to  pro¬ 
duce  seeds  (Elam  1994).  This  species  exhibits  low  seed  production,  which  is  attributed 
to  the  combined  effects  of  self-incompatibility  and  some  single-clone  populations  (Elam 
1994).  Likewise,  Indian  Knob  mountainbalm  is  expected  to  be  self-incompatible-.  Califor¬ 
nia  yerba  santa  is  a  fire-following  species.  Seeds  stored  in  the  soil  for  decades  germinate 
readily  during  the  first  spring  after  a  fire,  and  new  stems  sprout  from  underground  rhi¬ 
zomes  following  disturbances  such  as  fire  (Immel  2006).  Further,  in  the  absence  of  fire  to 
cue  seed  germination,  Eriodictyon  species  most  often  reproduce  vegetatively  by  rhizomes 
(USFWS  2013b).  Carlquist  et  al.  (1983)  described  species  in  the  genus  Eriodictyon  as  resis¬ 
tant  to  fire  and  drought  because  of  their  morphology  composed  of  underground  rhizomes 
and  aboveground  stems.  They  reported  Indian  Knob  mountainbalm  to  be  adapted  to  dry 
habitats,  with  the  underground  rhizomes  functioning  as  succulent  water-storage  and  stem 
innovation  organs. 

Wells  (1962)  suspected  that  the  rhizomes  of  Indian  Knob  mountainbalm  serve  for  re¬ 
production  when  the  aboveground  structures  are  destroyed  by  fire,  which  seems  likely. 
In  addition  based  upon  our  observations,  the  rhizomes  appear  to  be  the  primary  source 
of  reproduction  where  the  substrate  is  physically  disturbed,  such  as  by  a  rock  slide  or 


Knapp,  D.  2016.  2015  annual  report  to  U.S.  Fish  and  Wildlife  Service.  Santa  Barbara  Bot.  Garden, 
Calif..  16  pp. 
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Fig.  4.  Upper:  Ramets  (the  vertical  stems)  produced  from  underground  rhizomes  (indicated  by  black  ar¬ 
rows)  of  Indian  Knob  mountainbalm  Eriodictyou  altissiminn  at  southwest  Guidetti  Ranch,  San  Luis  Obispo 
County,  California  (occurrence  5),  5  May  2016.  Lower:  A  collapsed  individual  of  Indian  Knob  mountain¬ 
balm  growing  in  an  open  area  along  the  old  graded  dirt/rock  road  on  the  main  ridge  of  Indian  Knob 
mountain,  southwest  Guidetti  Ranch,  San  Luis  Obispo  County,  California  (occurrence  5),  5  May  2016. 
New  vertical  stems  have  sprouted  from  the  primary  stem,  which  is  now  horizontal  on  the  ground. 
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grading,  and  where  there  is  an  absence  of  recent  fire.  When  the  stem  is  removed  or  de¬ 
stroyed,  a  new  stem  may  be  produced  from  the  underground  rhizome  if  present.  We  ob¬ 
served  that  the  stems  of  Indian  Knob  mountainbalm  initially  stand  erect.  However,  because 
the  stems  and  branches  are  pliable,  taller  plants  bend  outward  when  not  supported  by  ad¬ 
jacent  vegetation  (Fig.  5).  This  results  in  some  stems  and  branches  touching  the  ground, 
and  eventually  collapsing  to  become  horizontal  on  the  ground  (Fig.  5).  This  presents  an 
ideal  situation  for  horizontal  stems  on  the  ground  to  be  covered  by  falling  rocks  or  shift¬ 
ing  sand  and  to  become  contact  points  for  layering  (subsequent  rooting  and  sprouting  of 
shoots).  We  observed  many  stems  in  this  horizontal  position  and  with  new  vertical  stems 
arising  from  the  horizontal  stems  (Fig.  4),  although  we  did  not  specifically  see  or  search 
for  new  roots.  Thus,  Indian  Knob  mountainbalm  appears  well-adapted  for  substrates  with 
physical  disturbances,  including:  steep  rocky  slopes  (Fig.  6),  cliff  faces,  fallen  rock  debris 
(Fig.  6),  sand  dunes  (shifting  sand;  Fig.  5),  roadsides  (Fig.  7),  old  graded  substrates  such 
as  dirt/rock  roads,  the  talus  of  graded  substrates,  and  trails. 

Indian  Knob  mountainbalm  can  reproduce  also  by  seeds,  which  are  minute  (~0.4  mm 
length,  mean  0.2  mg  mass;  Wells  1962).  Wells  (1962)  observed  a  large  production  of  seeds 
in  individuals  at  the  largest  occurrence  at  Indian  Knob  mountain.  However,  John  Chesnut 
(Los  Osos,  Calif,  pers.  comm,  to  Diane  Steeck,  USFWS,  1997)  observed  a  low  ratio  of 
seeds/ovules  (structures  that  develop  into  seeds  when  fertilized;  10%  or  less)  in  the  western 
occurrences  in  the  Los  Osos  area  and  also  a  low  number  of  seeds  per  ramet.  Some  occur¬ 
rences  may  have  low  genetic  diversity  due  to  rhizomatous  reproduction^  and  could  consist 
of  only  a  single  clone  (USFWS  1998,  2013a).  Although  the  role  of  fire  in  reproduction  is 
not  known,  it  obviously  reduces  dense  mature  vegetation  and  removes  debris  litter,  and  it 
may  be  needed  to  somehow  revitalize  mature  stands  of  Indian  Knob  mountainbalm^’'* 
(USFWS  1998)  or  to  cue  seed  production  or  germination  (USFWS  2013b).  USFWS 
(2013b)  discussed  the  potential  roles  of  fire  in  the  species’  ecology  in  a  status  report. 

As  of  2019,  Indian  Knob  mountainbalm  is  known  from  seven  occurrences,  with  a  total 
population  count  of  6,489+  individuals  in  2016  (Table  1).  The  seven  occurrences  comprise 
a  geographic  range  of  34  km-,  which  spans  a  distance  of  23  km  from  southeast  of  Indian 
Knob  (5. 1  km  south  of  the  city  of  San  Luis  Obispo),  northwest  to  Hazard  Canyon  (400  m 
south  of  Los  Osos)  in  Montana  de  Oro  State  Park  (Fig.  2).  This  area  has  a  Mediterranean 
climate  (warm  dry  summers,  cool  wet  winters).  The  dry  season  is  May  to  August,  with 
most  rain  falling  from  December  to  March.  Mean  annual  temperature  and  rainfall  are 
~15.1  and  ~536  mm,  respectively  (Ryan  1994).  We  document  the  occurrences  at  98 
to  263  m  elevation.  The  known  occurrences  are  at  the  two  ends  of  the  geographic  range, 
and  the  Irish  Hills  comprise  most  of  the  intervening  area.  We  suspect  the  species  occurs  at 
additional  locations  in  the  intervening  area,  however,  access  is  difficult  because  of  distance 
from  roads,  lack  of  trails,  steep  terrain  and/or  private  ownership. 

Six  occurrences  (1,  2,  3,  4,  6,  8)  of  Indian  Knob  mountainbalm  are  in  protected  areas 
(Table  2),  and  one  occurrence  (5)  is  mostly  in  a  protected  area.  Two  occurrences  (1,3)  are 
likely  extirpated,  one  occurrence  (5)  contains  a  large  number  of  plants  (6,346+  stems  in 
2016),  and  four  occurrences  (2,  4,  6,  8)  contain  low  numbers  of  plants  (^80  stems  in  2016) 
of  which  two  (occurrences  2,  4)  are  in  decline.  The  likely  extirpations  and  declines  may  be 
linked  to  the  altered  fire  regime,  specifically  the  lack  of  recent  fire  (fire  suppression)  in  the 


Bittman,  R.  1985.  Element  preservation  plan  for  Eriodictyon  altissitmnu  (Indian  Knob  mountain  balm). 
Nature  Conserv.,  San  Francisco,  Calif.,  5  pp. 
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Fig.  5.  Upper:  An  individual  of  Indian  Knob  mountainbalm  Eriodictyou  altissiinuin  growing  in  an  open 
area  on  the  main  ridge  of  Indian  Knob  mountain,  southwest  Guidetti  Ranch.  San  Luis  Obispo  County. 
California  (occurrence  5),  4  May  2016.  This  tall  plant  (3  to  4  m  height)  is  supported  behind  by  the  adjacent 
vegetation,  however,  it  is  collapsing  in  front  where  there  is  no  support.  Lower:  A  collapsed  individual  of 
Indian  Knob  mountainbalm  on  a  sand  dune  in  Morro  Dunes  Ecological  Reserve  East,  Los  Osos.  San 
Luis  Obispo  County,  California  (occurrence  6),  20  April  2016.  This  previously  tall  plant  (3  to  4  m)  is  now 
horizontal  on  the  ground  and  growing  in  an  open  area  along  the  trail.  Michael  J.  Walgren  of  California 
Department  of  Parks  and  Recreation  is  collecting  insects  while  David  Chipping  of  California  Polytechnic 
State  University  observes. 
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Fig.  6.  Upper:  A  colony  of  Indian  Knob  mountainbalm  Eriodiclyou  altissuinnu  growing  on  an  old 
rockslide  on  a  steep  slope  of  Indian  Knob  mountain,  southwest  Guidetti  Ranch.  San  Luis  Obispo  County, 
California  (occurrence  5),  5  May  2016.  Lower.  A  colony  of  Indian  Knob  mountainbalm  growing  in  fallen 
rock  debris  of  a  crumbling  cliff  face  on  Indian  Knob  mountain,  southwest  Guidetti  Ranch.  San  Luis  Obispo 
County.  California  (occurrence  5),  4  May  2016. 


Table  1 .  Reported  numbers  of  Indian  Knob  mountainbalm  Eriodictyon  altissimum  in  the  seven  known  occurrences  in  western  San  Luis  Obispo  County,  California. 
Occurrence  7  was  combined  into  occurrence  5  in  2018.  Occurrence  8  is  new.  X  =  species  observed.  Italics  indicate  the  words  as  used  by  the  footnoted  references. 
( Continued  on  next  page). 
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Table  2.  Conservation  status  of  the  seven  known  occurrences  of  Indian  Knob  mountainbalm  Eriodictyon  altissimum  in  western  San  Luis  Obispo  County,  California. 
Occurrences  8  is  new. 
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Fig.  7.  A  colony  of  Indian  Knob  mountainbalm  Eriodictyou  altissiniuni  beside  the  dirt/rock  road  in 
south  central  Guidetti  Ranch,  San  Luis  Obispo  County.  California  (new  occurrence  8),  1  December  2016. 
The  colony  is  1.14  km  southeast  of  the  nearest  colony  of  occurrence  5  on  southwest  Guidetti  Ranch  (below 
Indian  Knob).  The  stretch  of  road  at  occurrence  8  was  covered  with  soil/rock  from  below  Indian  Knob. 
It  seems  likely  the  plants  at  occurrence  8  originated  from  individuals  that  were  transported  with  soil/rock 
extracted  from  occurrence  5. 


habitat  (Gambs  and  Holland  1988;  USFWS  2013b).  Regarding  occurrence  5  that  is  mostly 
in  a  protected  area,  53%  (32  ha)  is  on  the  private  Guidetti  Ranch  with  a  conservation 
easement  to  the  City  of  San  Luis  Obispo,  40%  (24  ha)  is  on  multiple  private  properties 
in  Baron  Canyon  Ranch  Estates,  and  8%  (5  ha)  is  on  private  property  of  Pacific  Gas  and 
Electric  Company. 

The  primary  threat  to  the  species  identified  in  1994  has  been  eliminated:  habitat  loss  as 
a  result  of  potential  development  (USFWS  1994),  specifically  potential  surface  mining  at 
the  largest  occurrence  on  Indian  Knob  mountain,  and  potential  residential  development 
at  two  small  occurrences  in  Los  Osos.  At  the  private  Guidetti  Ranch  (occurrence  5; 
Tables  1,  2)  that  includes  much  of  Indian  Knob  mountain  and  where  the  bulk  of  the 
largest  population  exists,  the  City  of  San  Luis  Obispo  purchased  a  conservation  easement 
in  1 996  that  allows  continued  operation  of  the  ranch  for  livestock  grazing  while  protecting 
Indian  Knob  mountainbalm  and  its  habitat,  and  it  also  allows  limited  access  for  natural 
resource  conservation  and  public  education.  At  listing  in  1995,  occurrences  4  and  6  were 
on  two  private  properties  (83  ha)  in  Los  Osos  owned  by  a  development  corporation  that 
intended  to  build  houses.  These  two  properties  were  purchased  by  the  State  of  California 
in  2001,  and  they  now  comprise  CDFW’s  Morro  Dunes  Ecological  Reserve  East.  In  2019, 
dense  vegetation  is  a  threat  at  all  occurrences  except  occurrence  8  (see  details  for  each  oc¬ 
currence  in  Appendix  1).  Also  in  2019,  four  of  the  five  extant  occurrences  of  Indian  Knob 
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mountainbalm  are  threatened  by  stochastic  events  because  small  numbers  of  individuals 
comprise  them  (each  ^80  plants),  and  also  because  of  their  isolation  (distance)  from  other 
extant  occurrences.  These  threats  were  recognized  previously  by  USFWS  (1994). 

In  addition  to  the  immediate  threats  (Table  2),  climate  change  could  be  a  long-term 
threat  to  Indian  Knob  mountainbalm  because  California  is  becoming  hotter  and  drier. 
The  summers  of  2017  and  2016  were  each  the  warmest  in  California  since  record  keeping 
began  in  the  late  ISOO’s"^.  Considering  data  up  to  2015,  most  of  the  warming  occurred  in  the 
past  35  years  with  15  of  the  16  warmest  years  occurring  since  200C.  The  3-yr  period  from 
2012  to  2014  was  the  hottest  and  driest  in  California  in  the  100-yr  time  frame  considered  by 
Mann  and  Gleick  (2015),  and  it  was  the  most  severe  drought  in  California  in  the  past  1,200 
yrs  (Griffin  and  Anchukaitis  2014).  In  consideration  of  the  life  history  traits  used  to  assess 
climate  change  vulnerability  for  rare  plants  in  California  (Anacker  et  al.  2013),  Indian 
Knob  mountainbalm  would  be  considered  highly  vulnerable  to  climate  change.  These  traits 
include  low  dispersal  ability,  an  apparent  narrow  historical  climate  exposure,  dependence 
on  particular  disturbance  regimes,  habitat  specialization,  and  likely  low  genetic  diversity. 

Based  upon  our  field  observations,  the  several  statements  by  USFWS^^^  (2013a)  that 
invasive  nonnative  grasses,  particularly  perennial  veldt  grass  Ehrharta  calycina,  pose  the 
strongest  threat  to  Indian  Knob  mountainbalm  are  not  correct.  To  the  contrary,  no  oc¬ 
currence  of  Indian  Knob  mountainbalm  is  immediately  threatened  by  non-native  grasses. 
Although  perennial  veldt  grass  exists  in  the  surrounding  landscapes,  we  observed  it  only  at 
occurrence  5  with  just  a  few  clumps  at  one  location.  Although  perennial  veldt  grass  could 
become  a  threat  in  the  future  by  invading  an  occurrence  after  fire  (USFWS  2013a),  it  is 
more  abundant  in  coastal  sage  scrub  and  ruderal  and  disturbed  plant  communities  than  in 
chaparral. 

Using  our  information  and  international  standards  (lUCN  2012,  2014),  Indian  Knob 
mountainbalm  in  2019  meets  the  lUCN  criteria  for  endangered  including  the  following: 
geographic  range,  fragmented;  extent  of  occurrence,  34  km-  (<100  km-);  area  of  occu¬ 
pancy,  <2.3  km-  (<10  km^);  and  quality  of  habitat,  continuing  to  decline  (dense  vegeta¬ 
tion,  lack  of  recent  fire).  Considering  these  attributes,  Indian  Knob  mountainbalm  faces 
an  extremely  high  risk  of  extirpation  at  four  of  the  five  extant  occurrences.  Four  of  the 
five  extant  occurrences  are  small  with  ^80  individuals  each,  and  two  other  occurrences 
are  likely  extirpated.  One  occurrence  is  large  with  >6,000  individuals,  and  it  is  mostly  in  a 
protected  area. 


4 
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Recommendations. — Despite  increased  attention  in  recent  years,  little  is  known  about  the 
biology  and  ecology  of  the  endangered  Indian  Knob  mountainbalm  in  2019.  Therefore, 
coordinated  conservation  and  research  are  needed  to  further  understand  the  species,  and 
to  restore  and  maintain  the  five  extant  occurrences.  These  efforts  should  include  manage¬ 
ment  actions  to  benefit  the  occurrences,  searches  for  additional  locations,  studies  of  genetic 
diversity  and  reproductive  biology,  introduction  of  Indian  Knob  mountainbalm  into  liv¬ 
ing  collections  at  botanic  gardens,  investigations  of  potential  barriers  to  recruitment,  and 
investigations  of  the  species’  relationship  with  fire.  Specifically,  we  intend  to  make  a  con¬ 
servation  seed  collection  for  seed  banking  in  the  near  future,  and  a  subset  of  these  seeds 
will  be  used  to  investigate  germination  requirements,  including  cues  associated  with  fire¬ 
following  species.  In  addition,  cuttings  will  be  collected  and  cultivated  along  with  seeds  to 
develop  protocols  for  propagation  of  the  species.  Finally,  a  study  of  population  genetics  is 
planned  to  investigate  genetic  diversity  within  and  among  the  occurrences. 
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Appendix  1 

The  seven  known  occurrences  of  Indian  Knob  mountainbalm  Eriodictyon  altissimum  in  western  San  Luis 
Obispo  County  in  central  coastal  California. 

Occurrence  1.  ~35. 303567,  —120.846864;  119  m  elevation.  The  coordinates  that  we  give  are  in  a  former 
gully  ~1 19  m  E  of  the  intersection  of  Travis  Drive  and  Houston  Drive  in  Los  Osos.  The  location  is  on  the 
Broderson  site  (32  ha),  which  is  owned  by  the  County  of  San  Luis  Obispo  and  used  as  a  leach  field  (3  ha)  for 
recycled  water  from  the  Los  Osos  wastewater  treatment  plant.  The  occurrence  is  in  the  area  (29  ha)  protected 
by  deed  restriction  executed  in  2017,  and  managed  per  a  habitat  management  plan  approved  by  the  Cali¬ 
fornia  Coastal  Commission  in  2017  (Kate  Ballantyne,  Co.  of  San  Luis  Obispo,  Calif,  pers.  comm.  2017).  It 
is  on  Baywood  Fine  Sand  comprising  ancient  stabilized  sand  dunes  of  windblown  origin  during  the  late  to 
middle  Pleistocene  Epoch  (Carpenter  and  Storie  1928;  CDFG  2005;  Wiegers  2009).  We  searched  on  22  June 
2018,  10  April  2017  and  20  April  2016,  and  we  found  no  Indian  Knob  mountainbalm.  There  is  confusion 
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regarding  this  occurrence.  The  only  records  are  CDFG  (1979)  and  McLeod  (1986  in  CDFW  2017).  CDFG 
(1979)  states  the  following:  “Three  populations.  The  Indian  Knob  population  consists  of  12  individuals,  the 
Broderson  Road  site  includes  about  30,  and  the  Hazard  Canyon  population  is  slightly  larger.  Development 
activities  following  existing  high  volume  recreational  uses  have  caused  significant  decline  in  the  Roderson 
[sic]  Road  population.”  The  maps  of  CDFG  (1979)  and  McLeod  (1986  in  CDFW  2017)  both  depict  a  circle 
immediately  E  of  the  bend  of  Travis  Drive,  which  is  part  of  Cabrillo  Estates  (~250  homes).  McLeod  (1986 
in  CDFW  2017)  stated  the  following:  “Population  1  was  originally  discovered  by  Dirk  [Walters]  in  the  early 
’70s.  He  and  I  tramped  over  the  whole  area  last  year  and  couldn’t  locate  it.  Dirk  said  recently  that  it  might 
have  been  washed  away  as  it  was  along  a  gully.  Em  sorry  we  couldn’t  find  it  as  that  was  the  population  with 
which  you  would  probably  be  most  concerned.”  This  occurrence  and  occurrence  4  are  on  properties  that  in 
the  1980’s  were  proposed  for  the  community  sewage  system*^.  Gambs  and  Holland  (1988)  conducted  a  veg¬ 
etation  survey  of  the  properties  with  occurrences  1  and  4.  Although  Gambs  and  Holland  (1988)  discussed 
occurrence  4,  they  did  not  mention  occurrence  1.  Using  Google  Earth  aerial  imagery  (dated  27  May  1994), 
we  observed  that  a  gully  previously  existed  in  the  mapped  area  for  occurrence  1  (CDFG  1979;  McLeod  1986 
in  CDFW  2017)  and  parallel  to  Travis  Drive,  but  which  is  not  obvious  in  2017  because  of  dense  vegetation. 
Unconfirmed  word  of  mouth  is  that  the  Indian  Knob  mountainbalm  were  located  ~100  m  northeast  of  the 
intersection  of  Travis  Drive  and  Houston  Drive,  which  is  in  the  mapped  areas  of  CDFG  (1979)  and  McLeod 
(1986  in  CDFW  2017),  and  also  in  the  former  gully  reported  by  McLeod  (1986  in  CDFW  2017).  Because 
of  multiple  unsuccessful  attempts  to  locate  occurrence  1,  USFWS  (2013a,  b)  considered  occurrence  1  to  be 
extirpated.  In  contrast,  CDFW  considered  occurrence  1,  which  is  west  of  Broderson  Avenue  extended,  to 
be  a  mis-smapping  of  occurrence  4,  which  is  in  the  same  landscape  but  east  of  Broderson  Avenue  extended. 
Consequently,  CDFW  combined  occurrences  1  and  4  in  2013.  However,  instead  of  removing  occurrence  1 
from  the  California  Natural  Diversity  Database,  CDFW  removed  occurrence  4  while  retaining  occurrence 
1  but  with  coordinates  for  occurrence  4  (K.  Gross,  pers.  comm.  2017).  Although  no  herbarium  specimens 
document  occurrence  1,  it  appears  to  be  valid  as  indicated  by  CDFG  (1979)  and  McLeod  (1986  in  CDFW 
2017). 

Occurrence  2.  35.290381,  -120.859230;  35.290216,  -120.859599;  35.290156,  -120.859721;  35.291068, 

—  120.864935;  154  to  210  m  elevation.  The  coordinates  that  we  give  are  from  our  GPS  device  next  to  the 
plants.  Also,  Sayers  2008  (in  CDFW  2017;  Lazar  2018  pers.  comm.)  reported  two  additional  colonies  at 
the  following  two  locations  that  are  along  or  near  the  same  trail:  35.290202,  —120.860118;  35.290336, 

—  120.860566;  208  to  212  m.  We  surveyed  and  censused  this  occurrence  on  21  April  2016.  It  comprises 
six  disjunct  colonies  (we  refer  to  them  as  A,  B,  C,  D,  E  and  F,  from  east  to  west)  along  the  old  Manzanita 
Trail  on  a  south-facing  ridge  on  the  south  side  of  Hazard  Canyon  in  Montana  de  Oro  State  Park,  ~384  m 
south  of  Los  Osos.  The  trail  is  an  old  dirt/rock  road  between  the  east  section  of  the  Hazard  Peak  Trail  and 
the  Bloody  Nose  Trail  that  was  graded  likely  in  the  1940’s  or  1950’s  when  the  land  was  private  property. 
The  trail  is  closed  to  the  public.  The  location  of  this  occurrence  as  mapped  by  McLeod  (1981  in  CDFW 
2017)  is  off  by  >500  m.  The  greatest  number  of  recorded  individuals  for  this  occurrence  was  62  to  1 50  plants 
by  McLeod  (1985  in  CDFW  2017),  and  the  recorded  numbers  have  declined  since  then:  40  plants  in  1998 
(Hickson  and  Hillyard  1998  in  CDFW  2017),  37  plants  in  2005  (pers.  obs.),  28  plants  in  2008  (Sayers  2008 
in  CDFW  2017),  37  plants  in  2009  (pers.  obs.),  and  20  plants  in  2016  (pers.  obs.).  In  2016  this  occurrence 
comprised  20  plants  (15  big  individuals  with  woody  stems,  5  small  leafy  sprouts)  in  four  colonies  (A,  B,  C, 
F),  while  two  colonies  (D,  E)  reported  by  Sayers  (2008  in  CDFW  2017)  were  not  observed.  The  six  colonies 
are  separated  by  38  m,  10  m,  33  m.  44  m  and  402  m,  respectively,  spanning  a  straight-line  distance  of  530  m 
with  a  total  area  of  1 1,121  m-.  Colony  A  is  two  big  plants  at  the  top  of  a  cliff  face  (road  cut).  Sayers  (2008 
in  CDFW  2017),  likewise,  reported  two  plants  here  in  2008.  Colony  B  (six  large  stems,  one  little  stem)  and 
colony  C  (three  large  stems)  are  in  the  fallen  rock  debris  of  a  collapsed  cliff  face  (road  cut)  that  fell  in  201 1, 
and  the  bulk  of  the  two  colonies  were  covered.  Additional  plants  previously  grew  on  the  downslope  adja¬ 
cent  to  the  trail  and  were  covered  by  the  fallen  rock  debris.  The  rock  is  very  old  cemented  dune  sand  (D. 
Chipping,  pers.  comm.  2016),  and  it  is  friable  and  relatively  soft  (Wiegers  2009).  Colony  F  comprised  seven 
plants  (three  big  individuals  with  woody  stems,  four  small  leafy  sprouts)  in  an  area  of  38  m“  in  2016.  Four 
small  stems  were  growing  in  the  sandstone  road  (shale  according  to  McLeod  1985  in  CDFW  2017  is  incor¬ 
rect)  and  three  big  plants  in  the  rock  talus  pushed  down  the  slope.  Sayers  (2008  in  CDFW  2017)  reported 
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17  stems  here  in  2008,  and  McLeod  (1985  in  CDFW  2017)  reported  51-100  stems  in  1985.  With  only  seven 
stems  in  2016,  this  colony  is  in  decline  and  near  extirpation.  As  an  urgent  recovery  action  for  this  colony, 
we  recommend  that  the  vegetation  within  a  radius  of  8  m  be  thinned  and  that  the  potential  for  a  controlled 
miniburn  be  evaluated.  We  attribute  herbarium  specimen  CDA 17581,  which  was  collected  “ca.  2  mi  E  of 
Pecho  Road  along  S  leg  of  East  Hazard  Canyon  loop  trail.”  to  this  occurrence’^\  The  overland  distance 
from  Pecho  Valley  Road  to  occurrence  2  is  ~3.2  km  when  using  Hazard  Canyon  Road  and  the  ascending 
Manzanita  Trail. 

Occurrence  3.  ~35. 294023,  —120.850953;  239  m  elevation.  The  coordinates  that  we  give  are  from  our 
GPS  device  in  the  field.  This  occurrence  is  on  the  north  side  of  Hazard  Canyon  in  Montana  de  Oro  State 
Park  in  southwest  Los  Osos.  McLeod  (1981  and  1985  in  CDFW  2017)  reported  seven  plants  at  this  oc¬ 
currence  in  1981  and  1985,  but  no  one  has  reported  on  this  occurrence  since.  Although  McLeod  (1981  in 
CDFW  2017)  provided  a  map,  the  precise  location  has  remained  unknown.  We  searched  for  this  occurrence 
on  25  April  2017,  and  we  are  confident  that  we  found  the  general  location.  However,  we  saw  no  plants.  The 
soil  at  the  general  location  is  Baywood  Fine  Sand  on  a  ridge  top.  The  vegetation  is  now  dense,  predom¬ 
inantly  chamise  and  manzanita,  and  with  some  open  areas.  We  suspect  that  in  the  1980’s  the  vegetation 
was  less  dense,  and  in  the  ensuing  decades  the  vegetation  has  become  denser  and  crowded  out  the  Indian 
Knob  mountainbalm,  which  is  now  likely  extirpated  at  occurrence  3.  We  searched  also  along  Cable  Trail 
from  Pecho  Valley  Road  to  the  ridge  with  occurrence  3.  We  found  no  plants,  although  conditions  appeared 
suitable  in  several  places.  In  addition,  we  identified  a  nearby  area  (405  m  west  of  occurrence  3)  with  col¬ 
lapsed  clifi' faces  that  appeared  to  be  habitat.  We  tried  but  without  success  to  access  the  fallen  rock  debris 
from  Cable  Trail  above  and  from  Hazard  Canyon  Road  below.  Using  binoculars  we  saw  no  Indian  Knob 
mountainbalm.  In  brief,  unknown  colonies  of  Indian  Knob  mountainbalm  may  exist  along  the  north  side 
of  Hazard  Canyon. 

Occurrence  4.  35.305865,  -120.837545;  35.304063,  -120.838908;  35.302020,  -120.836210  (M.  Mooney, 
pers.  comm.  2018);  98  to  118  m  elevation.  The  coordinates  that  we  give  are  from  GPS  devices  next  to  the 
plants  in  three  colonies.  This  occurrence  is  in  southwest  Los  Osos  in  CDFW’s  Morro  Dunes  Ecological 
Reserve  East,  and  it  is  on  Baywood  Fine  Sand.  D.  Walters  identified  this  occurrence  in  1974  (McLeod 
1981  in  CDFW  2017).  Recorded  numbers  of  plants  were  six  individuals  in  1979  (McLeod  1981  in  CDFW 
2017),  30  in  1985  (McLeod  1985  in  CDFW  2017),  25  in  1988  (Gambs  and  Holland  1988),  and  23  in  2016 
(pers.  obs.).  There  is  confusion  regarding  this  occurrence,  whose  specific  location  became  “lost”  after  1988. 
McLeod  (1981  in  CDFW  2017)  mapped  the  location  of  occurrence  4  at  73  to  85  m  elevation  between 
the  sand  extensions  of  Palisades  Avenue  and  Ravenna  Avenue,  but  closer  to  the  former.  McLeod  (1985  in 
CDFW  2017)  described  the  location  as  “south  of  the  end  of  Palisade  St  Los  Osos  100  yds  W  of  stand  of 
Bishop  pine  [Finns  miiricatciY  at  91  m  elevation  but  marking  79  to  85  m  elevation  on  his  map.  Gambs  and 
Holland  (1988)  mapped  the  location  of  Bishop  pine  in  the  landscape,  showing  two  stands  near  Palisades 
Avenue  sand  extension  and  the  most  relevant  one  at  98  m  elevation  (35.304162,  —120.837298).  Gambs  and 
Holland  (1988)  described  the  location  of  occurrence  4  as  two  groups:  one  group  of  25  individuals  ~100  m 
west  of  the  sand  extension  of  Palisades  Avenue  at  ~76  m  elevation,  and  “an  isolated  individual  discovered 
by  Dirk  Walters  in  the  early  1970’s  south  of  the  end  of  Ravenna  Avenue”  that  they  could  not  find.  CDFG 
(2006)  recognized  occurrence  4  with  coordinates  35.30563,  —120.83734  at  83  m  elevation.  We  are  not  able 
to  reconcile  the  previously  stated  elevations  and  maps  for  occurrence  4  with  our  own  observations.  Because 
of  the  multiple  unsuccessful  attempts  to  find  this  occurrence,  USFWS  (2013a,  b)  considered  occurrence 
4  to  be  extirpated.  Also  in  2013,  CDFW  combined  occurrence  4  with  occurrence  1  because  they  believed 
that  occurrence  1  (35.30467,  —120.84598;  CDFG  2006)  was  a  mapping  error.  CDFW  removed  occurrence 
4  from  the  California  Natural  Diversity  Database  and  retained  occurrence  1  but  with  the  new  coordinates 
35.30556,  —120.83771  at  79  m  elevation,  which  is  35  m  southwest  of  the  former  occurrence  4  (CDFG  2006). 
However,  we  affirm  that  our  occurrence  4  is  separate  and  distinct  from  our  occurrence  1 .  J.  Chesnut  of  the 
California  Native  Plant  Society  (San  Luis  Obispo)  led  us  to  occurrence  4  on  20  April  2016.  We  observed 
23  individuals  in  two  colonies  separated  by  40  m:  colony  A  with  21  plants  (10  big  individuals  with  woody 
stems,  1 1  small  sprouts;  14  m“)  in  a  clearing,  and  a  more  western  colony  B  with  two  big  individuals  with 
woody  stems  (3  m“)  at  the  edge  of  California  live  oak  woodland  Querciis  agrifolia  with  ceanothus  Ceanothus, 
manzanita  Manzanita  and  chamise  Adenostonui  fasciculatiun.  J.  Chesnut  (pers.  comm.  2016)  had  observed 
the  two  colonies  previously  during  field  work  in  1986,  and  he  recalled  there  were  other  individuals  on  the 


Consortium  of  California  Herbaria.  2017.  UC/JEPS:  Consortium  search  results  for  Eriodictyon  al- 
tissinnun.  Available  (Nov.  2017):  http://ucjeps.berkeley.edu/consortium,  2  p. 
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sand  ridge  extending  north.  We  subsequently  searched  the  sand  ridge  on  10  April  2017  but  saw  no  other 
plants.  We  observed  colony  B  again  on  10  April  2017,  and  one  plant  was  dead  and  the  other  dying.  The 
coordinates  that  we  give  for  colony  A  are  139  m  west  of  the  stand  of  Bishop  pine  at  98  m  elevation  near 
the  end  of  the  sand  extension  of  Palisades  Avenue  (Gambs  and  Holland  1988),  which  is  likely  the  location 
referred  to  by  McLeod  (1981  and  1985  in  CDFW  2017)  and  Gambs  and  Holland  (1988).  Colony  B  is  at  the 
edge  of  dense  vegetation,  and  it  is  40  m  west  of  colony  A  and  even  closer  to  the  sand  extension  of  Ravenna 
Avenue.  In  August  2018,  M.  Mooney  (pers.  comm.  2018)  observed  one  individual  ~30  cm  height  at  a  new 
location  (colony  C),  330  m  southeast  of  colony  A  and  in  chaparral.  Five  herbarium  specimens  document 
occurrence  4:  CHSC69826and  HSC77146  collected  in  1982;  and  RSA522505,  SBBG95658  and  UC1583847 
collected  in  1986'*\  The  total  area  is  0.62  ha. 

Occurrence  5.  Based  upon  the  information  for  occurrences  5  and  7  in  a  draft  of  this  paper  that  was 
reviewed  by  CDFW,  in  2018  they  officially  combined  occurrence  7  into  occurrence  5  because  the  two  oc¬ 
currences  were  no  longer  separated  by  >0.4  km  (K.  Lazar,  pers.  comm.  2018).  However,  because  of  the 
different  land  ownerships  and  different  management,  we  present  the  separate  information  for  occurrences 
5  and  7  prior  to  2018  (see  also  Table  1).  The  total  known  occupied  area  is  ~60  ha:  Guidetti  Ranch,  32  ha; 
Baron  Canyon  Ranch  Estates.  24  ha;  and  Pacific  Gas  and  Electric  Company,  5  ha.  The  approximate  center 
of  occurrence  5  is  now  35.200234,  —120.669741.  The  total  count  in  2016  was  6,346+  stems.  USFWS  (1994) 
reported  350  stems  in  1991. 

Occurrence  5  before  2018.  35.200538,  -120.669342  to  35.198472,  -120.657299;  and  35.201136, 

—  120.662712  to  35.196987,  —120.662644;  197  to  263  m  elevation.  The  coordinates  that  we  give  are  from 
our  field  observations  and  Google  Earth  aerial  imagery.  The  approximate  center  of  the  occurrence  is  at 
35.199150,  —120.663379.  This  occurrence  is  narrowly  situated  on  Indian  Knob  mountain  and  the  nearby 
ridges  and  upper  slopes,  primarily  on  the  private  Guidetti  Ranch  and  extending  south  onto  property  of 
Pacific  Gas  and  Electric  Company.  The  substrate  is  sandstone  (Chipping,  pers.  comm.  2016).  We  surveyed 
and  censused  this  occurrence  on  4  and  5  May  2016.  We  counted  5.720  stems  (4,610  woody  stems,  1,110 
leafy  sprouts;  and  205  dead  stems)  on  the  Guidetti  Ranch  and  adjacent  property  of  Pacific  Gas  and  Electric 
Company.  This  is  a  minimum  number  because  there  were  many  areas  we  did  not  search  due  to  steepness 
of  slopes  or  lack  of  paths  through  dense  vegetation.  The  plants  are  primarily  on  the  ridges  and  south¬ 
facing  slopes.  Previous  reports  regarding  numbers  of  plants  were  12  in  1979  (CDFG  1979),  >100  in  1985 
(McLeod  1985  in  CDFW  2017),  and  200  in  2015“.  Our  substantially  higher  number  is  likely  the  result  of 
greater  effort,  using  a  systematic  method,  and  the  counting  of  every  stem.  The  known  occupied  area  com¬ 
prises  at  least  32  ha,  which  is  mostly  on  the  Guidetti  Ranch  (30  ha)  but  also  extending  south  and  west 
onto  private  property  owned  by  the  Pacific  Gas  and  Electric  Company  (3  ha).  Wells  (1962)  described  Indian 
Knob  mountainbalm  on  Guidetti  Ranch  as  a  pioneer  species  in  chaparral  dominated  by  Santa  Margarita 
manzanita  Arctostaphylos pilosula,  being  aggressive  along  roads  and  confined  to  shallow  sandy  soils  on  the 
sandstone  ridges.  CDFW  (2017)  described  the  microhabitat  as  sandstone  ridges  in  open  disturbed  areas  in 
chaparral  at  90  to  270  m  elevation.  In  addition  to  the  roadsides,  old  graded  dirt/rock  roads  and  ridges,  we 
observed  the  species  in  the  open  areas  on  the  steep  slopes  usually  with  a  disturbed  substrate,  and  also  in 
the  talus  of  graded  dirt/rock  roads,  on  cliff  faces,  in  the  fallen  rock  debris  of  collapsed  cliff  faces,  on  rock 
slides,  and  along  trails.  We  searched  several  disturbed  areas  in  the  chaparral  created  by  tree  falls,  which 
were  characterized  by  a  thin  layer  ( 1 5  cm)  of  fine  organic  debris  over  sandy  soil,  but  found  no  Indian  Knob 
mountainbalm  there.  On  Guidetti  Ranch  it  appears  that  Indian  Knob  mountainbalm  is  adapted  for  areas 
with  disturbed  soil  and  rock  where,  in  the  absence  of  recent  fire,  reproduction  seems  most  likely  by  rhizomes. 
We  suspect  that  reproduction  by  seeds  may  need  fire,  and  there  have  been  no  wildfires  on  Guidetti  Ranch 
since  before  the  1950’s  (Terry  Guidetti,  pers.  comm.  2016). 

Occurrence  7  before  2018.  35.209470,  -120.677157  to  35.200040,  -120.670969;  and  35.203853, 

—  120.670995  to  35.205170,  —120.676575;  193  to  235  ni  elevation.  The  coordinates  that  we  give  are  based  on 
our  field  observations  and  using  Google  Earth  aerial  imagery.  The  known  occupied  area  comprises  ~28  ha, 
with  the  center  at  35.203840,  —120.674061 .  This  occurrence  is  mostly  on  Baron  Canyon  Ranch  Estates  (277 
ha),  which  is  a  private  gated  community  in  the  mountains  and  hills  2.25  km  south  of  the  city  of  San  Luis 
Obispo.  The  community  comprises  29  parcels  of  land  for  houses  (0.4  to  34.0  ha)  and  three  common  areas 
(172  ha).  Most  houses  are  on  the  ridges  and  hills,  which  are  mostly  undeveloped  and  with  predominantly 
native  vegetation.  The  section  of  the  primary  mountain  ridge  with  Indian  Knob  mountainbalm  is  oriented 
mostly  in  a  northwest  to  southeast  direction,  and  Balm  Ridge  Way  is  the  paved  road  that  extends  along 
it.  This  occurrence  of  Indian  mountainbalm  is  on  the  primary  mountain  ridge  and  adjacent  slopes,  span¬ 
ning  a  distance  of  1.2  km.  It  is  on  the  same  primary  mountain  ridge  as  occurrence  5,  which  is  mostly  on 
Guidetti  Ranch  and  immediately  southeast  of  Baron  Canyon  Ranch  Estates.  Occurrence  7  extends  for  a 
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short  distance  S  onto  property  owned  by  Pacific  Gas  and  Electric  Company  and  east  onto  Guidetti  Ranch. 
The  southern  part  of  occurrence  7  adjoins  occurrence  5.  We  observed  Indian  Knob  mountainbalm  primar¬ 
ily  in  a  belt  of  Santa  Margarita  manzanita  that  appeared  to  stretch  continuously  across  the  several  ridges 
and  slopes.  We  counted  626  plants  (598  woody  stems,  28  leafy  sprouts)  on  27  July  2016,  and  there  were 
likely  many  more  plants  that  we  did  not  observe  because  of  the  large  occupied  area,  the  steepness  of  some 
slopes,  and  the  lack  of  access  through  the  dense  vegetation.  We  saw  many  tall  individuals  of  Indian  Knob 
mountainbalm  growing  within  and  rising  above  the  dense  mature  chaparral,  including  one  individual  that 
was  ~5.5  m  tall.  Similarly,  Wells  (1962)  observed  Indian  Knob  mountainbalm  “often  overtopping  by  five 
feet  or  more  the  even-statured  young  manzanitas  dating  from  the  last  chaparral  fire.”  LynneDee  [Oyler]  Al- 
thouse  (1991  in  CDFW  2017;  Paso  Robles,  Calif,  pers.  comm.  2017)  reported  >500  plants  (350  flowering, 
150  vegetative)  in  1991.  Current  threats  include  clearing  of  the  native  vegetation  by  private  landowners.  In 
2016  we  observed  illegal  clearing  of  at  least  310  m“  of  native  vegetation  on  one  private  property  contrary 
to  ordinance  of  San  Luis  Obispo  County,  and  it  included  destruction  of  some  Indian  Knob  mountainbalm, 
which  was  being  investigated  by  CDFW. 

Occurrence  6.  35.301508,  —120.830263;  35.301639,  —120.830129;  115  to  118  m  elevation.  This  occur¬ 
rence  is  in  south  Los  Osos  in  Morro  Dunes  Ecological  Reserve  East  on  Baywood  Fine  Sand.  It  is  25  m 
northwest  of  the  white  water  tank  near  Calle  Cordoniz  Road.  The  coordinates  are  from  our  GPS  device 
next  to  the  plants.  We  surveyed  and  censused  this  occurrence  on  20  April  2016,  and  we  counted  20  plants: 
1 1  big  individuals  with  woody  stems,  and  nine  small  sprouts.  The  plants  are  in  two  colonies  (20  m“,  32  m-) 
separated  by  1 1  m,  and  with  a  total  area  of  63  m-.  Much  of  Los  Osos  is  underlain  by  old,  stabilized  sand 
dunes  (Wiegers  2009),  and  this  occurrence  is  near  the  crest  of  a  sand  dune.  The  plants  are  along  several 
well-established  trails  and  at  the  edge  of  dense  mature  chaparral.  Several  of  the  taller  individuals  had  col¬ 
lapsed  and  were  growing  with  the  primary  stems  horizontal  across  the  ground.  Previous  records  for  this 
occurrence  are  15  plants  in  2012  (USFWS  2013b),  10  in  2010  (Butterworth  2010  in  CDFW  2017),  20  to  25 
in  2009  (USFWS  2013b),  and  1 1  to  50  plants  in  1985  (McLeod  1985  in  CDFW  2017).  We  identify  two  pri¬ 
mary  threats  to  this  occurrence:  the  dense  vegetation  immediately  surrounding  the  plants,  and  the  removal 
of  Indian  Knob  mountainbalm  by  hikers  and  equestrians  who  trim  the  vegetation  in  an  effort  to  keep  the 
trails  passable  (pers.  obs.  2017).  Therefore,  we  recommend  that  CDFW  (or  their  authorized  agents)  thin 
the  vegetation  in  vicinity  of  the  Indian  Knob  mountainbalm,  and  also  place  educational  signs  at  strategic 
locations  to  inform  hikers  and  equestrians  about  sensitive  plants  in  the  habitat. 

Occurrence  8  (new).  35.190008,  —120.650146;  124  m  elevation.  The  coordinates  that  we  give  are  from 
Google  Earth  aerial  imagery.  While  working  with  vegetation  contractors  of  the  Pacific  Gas  and  Electric 
Company  on  1  December  2016,  we  encountered  a  disjunct  colony  of  Indian  Knob  mountainbalm  (78  woody 
stems,  2  leafy  sprouts)  that  is  1.14  km  southeast  of  the  nearest  colony  of  occurrence  5  below  Indian  Knob  on 
Guidetti  Ranch.  The  new  occurrence  is  beside  the  dirt/rock  road  and  in  the  south  central  part  of  Guidetti 
Ranch  (Fig.  7),  which  is  in  the  area  protected  by  conservation  easement  to  the  City  of  San  Luis  Obispo.  The 
soil/rock  at  occurrence  8  and  occurrence  5  below  Indian  Knob  appear  similar  using  Google  Earth  aerial 
imagery.  Although  the  landscape  between  the  two  occurrences  is  intact  and  mostly  vegetated,  there  appears 
to  be  no  similar  soil/rock  in  this  intervening  area.  The  location  below  Indian  Knob  was  previously  surface 
mined,  and  the  extracted  soil/rock  was  used  locally  for  road  fill  and  road  repair.  T.  Guidetti  (pers.  comm. 
2018)  stated  that  the  stretch  of  road  at  occurrence  8  was  covered  with  soil/rock  from  below  Indian  Knob. 
Thus,  it  seems  likely  the  plants  at  occurrence  8  originated  from  individuals  that  were  transported  with 
soil/rock  extracted  from  occurrence  5.  Threats  to  this  occurrence  include  low  number  of  plants,  vehicles 
running  over  the  plants,  and  road  maintenance. 
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Optimizing  a  Municipal  Wastewater-based  Chlorella  vulgaris 
Photobioreactor  for  Sequestering  Atmospheric  CO2 
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Abstract. — Microalgae  photobioreactors  are  among  the  most  effective  systems  for 
capturing  gaseous  CO2,  the  main  contributor  to  global  warming.  Their  eapacity  to 
generate  massive  amounts  of  biomass  has  been  exploited  serendipitously  to  sequester 
CO2  and  explieitly  to  remove  nutrients  from  municipal  wastewater.  Unfortunately,  re¬ 
search  in  this  area  has  not  included  merging  these  dual  eapacities  to  address  global 
warming.  Instead,  most  are  focused  on  thermolytie  eonversion  of  biomass  into  energy 
whieh  in  end  returns  CO2  to  the  atmosphere.  In  this  study,  we  investigated  the  poten¬ 
tial  of  combining  the  two  microalgal  capacities  (that  of  deriving  nutrients  from  munie- 
ipal  wastewater  and  metabolic  carbon  from  atmospheric  or  industrial  CO2  supplies), 
into  an  integrated  means  of  redueing  nutrients  in  oeean-bound  wastewater  and  CO2 
in  the  atmosphere  simultaneously.  The  test  species  used  in  this  study  was  Chlorella 
vulgaris  (C.  vulgaris)',  the  turbidity  of  C.  vulgaris  was  used  as  a  measure  of  yield  in 
biomass.  Our  results  show  (i)  that  an  open  photobioreaetor,  and  not  a  closed  one,  is 
the  most  productive,  espeeially  when  augmented  with  industrial  CO2  (hence  making 
a  strong  ease  for  scrubbing  CO2  gas  from  industrial  sources),  (ii)  that  a  meehanieally 
agitated  C.  vulgaris  culture  is  more  productive  than  a  static  one,  (iii)  that  without  me¬ 
chanical  agitation,  32  ib  3  days  of  incubation  are  needed  to  reach  the  maximum  yield 
of  an  open  photobioreaetor,  (iv)  that  the  optimal  proportion  of  wastewater  (%WW) 
required  to  support  C.  vulgaris  growth  is  80  ±  3%;  at  least  33%  WW  is  required  to  ob¬ 
serve  growth  above  background,  and  (v)  that  without  intervention,  the  upper  pH  limit 
of  a  WW-based  C.  vulgaris  culture  is  8.69  ±  0.09.  Two  mutually  independent  models 
are  proposed  to  aide  in  sealing  up  an  open  WW-based  C.  vulgaris  photobioreaetor. 


Global  warming  and  eutrophication  are  two  of  the  most  pressing  challenges  confronting 
environmental  scientists  in  recent  times  (Committee  on  Science  for  EPA’s  Future  2012). 
Global  warming,  the  steady  heating  of  the  earth  surface  and  of  the  atmosphere,  is  partly  a 
consequence  of  atmospheric  emission  of  greenhouse  gases  such  as  carbon  dioxide  (CO2), 
methane  (CH4),  nitrous  oxide  (N2O),  and  some  fluorinated  compounds  (Rey  et  al.  2018). 
These  greenhouse  gases  do  not  only  absorb  and  retain  infra-red  energy  (Rey  et  al.  2018), 
they  prevent  the  release  of  heat  into  the  outer  space  from  the  earth  surface  (Akitt  2018). 
Of  these  gases,  CO2  is  implicated  the  most  in  global  warming  (reviewed  by  Anderson  et  al. 
2016).  In  2015  alone,  82%  of  all  greenhouse  gases  emitted  was  C02.^  In  modern  times. 


*  Corresponding  author:  ochan.otim@lacity.org.  ORCID:  https://orcid.org/0000-0001-7272-4356 
'  Olivier  J.G.J.,  G.  Janssens-Maenhoiit,  M.  Muntean  and  J.A.H.W  Peters.  2016.  Trends  in 
global  CO2  emissions.  The  Hague:  PBL  Netherlands  Environmental  Assessment  Agency.  Available 


42 


WASTEWATER-BASED  CHLORELLA  VULGARIS  BIOREACTOR  FOR  SEQUESTERING  CO.  43 


the  main  source  of  CO2  is  the  combustion  of  fossil  fuel  to  meet  domestic  and  industrial 
energy  needs  (Fischer  et  al.  2018;  Carotenuto  et  al,  2018).  Efforts  to  curb  global  warming 
henceforth  will,  therefore,  have  to  include  controlling  this  anthropogenic  CO2  emission 
alongside  reducing  CO2  levels  already  present  in  the  atmosphere. 

Eutrophication,  the  excessive  plant  growth  in  large  water  bodies,  is  due  to  increased 
availability  of  sunlight,  carbon  dioxide,  and  nutrients  (Greeson  1969;  Schindler  2006).  This 
excessive  plant  growth  in  turn  creates  areas  in  aquatic  environment  where  the  amount  of 
dissolved  oxygen  is  not  enough  to  support  marine  life.  In  current  times,  the  discharge  of 
large  amounts  of  nitrogen  and  phosphorus  into  the  aquatic  ecosystems  has  increased  the 
rate  of  eutrophication  and  the  extent  to  which  eutrophication  is  affecting  inland  waterways 
and  coastal  waters  (Carpenter  et  al.  1998;  Chislock  et  al.  2013).  The  affected  areas  (also 
referred  to  as  dead  zones)  have  been  shown  to  occur  frequently  near  where  agricultural 
runoff  and  wastewater  treatment  effluent  are  released  into  large  water  bodies.-  Efforts  to 
remediate  eutrophication  resulting  from  these  anthropogenic  activities  may  have  to  include 
polishing  agricultural  runoffs  and  municipal  wastewater  before  discharge. 

Over  the  last  few  years,  microalgae  have  emerged  as  promising  vehicles  for  converting 
CO2  into  large  quantity  of  biomass.  Research  in  this  area  though  is  primarily  focused  on 
optimizing  the  biomass  production  for  the  biofuel  industry  with  little  considerations  be¬ 
ing  given  to  the  consequential  release  of  the  CO2  back  into  the  atmosphere  (Xu  et  al. 
2006;  Sayre  2010;  Mondal  et  al.  2017).  In  our  opinion,  the  biomass  so  produced  could  be 
stored  to  control  the  amount  of  gaseous  CO2  in  circulation.  Storage  could  mean  delay¬ 
ing  or  preventing  the  release  of  CO2  into  the  atmosphere.  Others  have  proposed  chemical 
means  (Obersteiner  et  al.  2001;  Lackner  et  al.  2012;  Wurzbacher  et  al.  2012)  and  the  ocean 
(Shaffer  2010)  as  solutions  towards  reducing  the  amount  of  CO2  in  circulation.  However, 
chemical  means  would  require  careful  attention  to  the  disposal  of  the  chemical  wastes  so 
generated  in  the  process,  and  both  legal  and  environmental  hurdles  will  have  to  be  over¬ 
come  before  anthropogenic  CO2  is  stored  in  the  ocean. 

In  this  study,  we  propose  a  natural  way  of  sequestering  CO2  in  biomass.  We  investigate 
whether  a  stable  culture  of  microalgae  with  the  following  characteristics  can  be  maintained 
in  a  mixotropic  photobioreactor:  (i)  one  that  derives  metabolic  carbon  from  CO2  (Sayre 
2010;  Mondal  et  al.  2017;  Raven  2017)  and  (ii)  one  that  would  also  utilize  secondary- 
treated  municipal  wastewater  as  the  only  source  of  nutrients  (Lau  et  al  1995;  Wang  et  al. 
2010;  Craggs  et  al.  2011;  Shi  et  al.  2016).  We  believe  that  a  successful  integration  of  these 
two  microalgal  natural  capacities  can  play  a  role  in  scrubbing  CO2  gas  not  only  from  the 
atmosphere  but  from  industrial  sources  as  well  while  removing  nutrients  from  ocean-bound 
municipal  wastewater  (WW)  at  the  same  time.  This,  metaphorically,  is  killing  two  birds  with 
the  same  stone.  The  removal  of  nutrients  in  particular  is  of  importance  to  watersheds  and 
municipal  wastewater  treatment  plants  that  discharge  freshwater  and  treated  wastewater 
effluents,  respectively  into  large  water  bodies.  Elevated  level  of  anthropogenic  nutrients  in 
the  ocean  for  example  is  known  to  promote  toxic  algal  blooms  (Anderson  et  al.  2002). 

The  key  factors  investigated  in  this  study  are  those  known  to  affect  the  productivity 
of  microalgae  photobioreactor.  These  factors  include  microalgae  species  to  culture,  the 


from:  http:/ / edgar.jrc.ec.europa.eu/news_docs/jrc-2016-trends-in-global-co2-emissions-2016-report- 
103425.pdf.  Accessed  05  July,  2018. 

Paine  V.  2012.  What  Causes  Ocean  “Dead  Zones”?  Sci.  Am.,  Available  from:  https://www. 
scientificamerican.com/article/ocean-dead-zones/.  Accessed  15  January,  2019. 
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Table  1.  Composition  of  secondary-treated  wastewater  used  in  this  study.  ^ 


Test 

Results 

pH 

7.3 

Carbonaceous  Biochemical  Oxygen  Demand 

5  mg/L 

Biological  Oxygen  Demand 

12mg/L 

Settleable  Solids 

<0.1  mL/L 

Turbidity 

5.1  NTU 

Mercury,  Hg 

0.00239  |ig/L 

Silver.  Ag 

0.061  |ig/L 

Boron,  B 

0.541  mg/L 

Antimony,  Sb 

2.75  pg/L 

Beryllium,  Be  ‘‘ 

0.03  |ig/L 

Cadmium,  Cd  ‘‘ 

0.05  pg/L 

Copper,  Cu 

13.4  pg/L 

Nickel,  Ni 

6.52  pg/L 

Thallium,  Th  ‘‘ 

0.02  pg/L 

Zinc.  Zn 

14.8  pg/L 

Arsenic,  As 

2.30  pg/L 

Selenium,  Se  ^ 

0.52  pg/L 

Total  Suspended  Solids 

9.7  mg/L 

Total  Dissolved  Solids 

768  mg/L 

Total  Organic  Carbon 

19.1  mg/L 

Total  Kjeldahl  Nitrogen 

47.8  mg/L 

Organic  Nitrogen 

<1  mg/L 

Total  Phosphate 

2.82  mg/L 

Ammonia  as  Nitrogen 

47.0  mg/L 

Nitrate  as  Nitrogen 

<0.1  mg/L 

Nitrite  as  Nitrogen  ‘‘ 

0.32  mg/L 

Chloride 

303  mg/L 

^  Available  data  at:  http://ciwqs.waterboards.ca.gov/ciwqs/readOnly/PublicReportEsmrAtGlance 
Servlet?reportID  =  2&isDrilldown  =  true&documentID  =  1893752. 

Values  lie  between  detection  limit  and  regulatory  reporting  limit. 


method  of  delivering  CO2  to  the  culture,  the  nutritional  needs  of  microalgae,  the  pH 
changes  occurring  during  incubation  in  the  microalgae  culture,  and  photobioreactor  con¬ 
figuration  (i.e.,  an  open-air  system  versus  a  closed  one;  most  researchers  use  the  latter 
without  provision  for  releasing  metabolic  O2  known  to  be  toxic  to  microalgae).  Earlier 
studies  in  our  lab  (unpublished)  suggests  that  C.  vulgaris  can  flourish  under  the  mixotropic 
cultivation  setup  adopted  for  this  study.  The  potential  of  using  C.  vulgaris  similarly  has 
been  demonstrated  recently  by  Shi  et  al.  (2016),  and  the  kinetics  of  C  vulgaris  growth  in 
freshwater  has  also  published  by  Adamczyk  et  al.  (2016). 

Materials  and  Methods 

All  containers  used  in  this  study  were  thoroughly  washed  with  soapy  water,  rinsed  with 
distilled  water  (DH2O)  and  disinfected  with  3%  aqueous  solution  of  hydrogen  peroxide 
(PL  Development,  Clinton,  SC,  USA).  Secondary-treated  wastewater  was  acquired  from 
the  Hyperion  Water  Reclamation  Plant  (City  of  Los  Angeles,  CA,  USA);  its  parameters, 
listed  in  Table  1,  were  determined  as  described  elsewhere  (Otim  et  al.  2018). 

Chlorella  vulgaris  cultures  (10  mL  in  Alga-Gro  Lreshwater  growth  medium)  were  pur¬ 
chased  from  Carolina  Biological  Supply  (Burlington,  NC,  USA)  and  cultured  using  a 
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Fig.  1.  (a)  The  setup  used  to  maintain  a  viable  stock  of  Chlorella  vulgaris  in  distilled  water,  (b)  Cul¬ 

tivating  Clilorella  vulgaris  in  secondary-treated  wastewater  (WW)  collected  from  the  City  of  Los  An¬ 
geles  Hyperion’s  5-Mile  Outfall  system.  The  wastewater  was  serially  mixed  (in  inL  of  WW)  with  dis¬ 
tilled  water  to  a  final  volume  of  310  inL  in  each  500-mL  flask  (or  500-mL  beaker,  not  shown)  as  fol¬ 
lows:  Flask  1:  0  (0%  WW),  Flask  2:  50  (16%  WW),  Flask  3:  100  (33%  WW),  Flask  4:  150  (50%  WW), 
Flask  5:  200  (67%  WW),  Flask  6:  250  (83%  WW),  aud  Flask  7:  300  (100%  WW).  The  310  mL  total  vol¬ 
ume  includes  10  inL  of  laboratory  stock  of  Chlorella  vulgaris  (see  Materials  and  Methods  section).  Glass 
wool  and  perforated  saran  wraps  were  used  to  cover  the  open  mouths  of  Erlenmeyer  flasks  and  beakers, 
respectively. 


modification  of  James  (2012)  protocol  as  follows.  The  C  vulgaris  cultures  were  first  stored 
as  received  at  room  temperature  under  an  alternating  cycle  of  artificial  white  light  (16  h) 
and  darkness  (8  h)  for  seven  days  to  allow  acclimatization.  Then  before  use,  the  cultures 
were  diluted  to  1  L  with  fresh  Alga-Gro  Freshwater  growth  medium  and  incubated  for 
seven  more  days  under  alternating  cycle  of  artificial  white  light  as  described  above.  To 
maintain  a  laboratory  stock,  C.  vulgaris  was  cultured  in  an  open  photobioreactor  (30- 
L  DH2O)  fitted  with  the  capacity  to  boost  the  level  of  CO2  in  microalgae  culture  au¬ 
tomatically  via  a  gas  dispersion  frit  to  bring  down  pH  to  7.0  when  necessary  (Fig.  la); 
the  upper  limit  of  pH  was  set  to  7.4.  Nutrients  (20  mL  of  Guillard’s  (F/2)  Marine  Wa¬ 
ter  Enrichment  Solution;  Carolina  Biological  Supply,  Burlington,  NC;  Guillard  and  Ry- 
ther  1962;  Guillard  1975)  were  supplied  to  the  C.  vulgaris  laboratory  stock  every  seven 
days. 

To  determine  the  ratio  of  wastewater-to-distilled  water  (WW/DH2O)  suitable  for  opti¬ 
mal  growth,  two  separate  trials  of  seven  310-mL  photobioreactors  were  set  up  in  500-mL 
beakers  (not  shown)  and  in  Erlenmeyer  flasks  shown  in  Fig.  lb  (labeled  1  to  7)  to  as¬ 
sure  reproducibility.  This  dual  setup  also  served  to  gauge  in  a  rudimentary  way  the  effect 
of  exposing  different  surface  areas  on  photobioreactor  performance.  The  WW/DH2O  ra¬ 
tios  were  formulated  as  follows  (v/v):  0/300,  50/250,  100/200,  150/150,  200/100,  250/50, 
300/0.  Into  each  photobioreactor  was  added  10  mL  of  laboratory  stock  C.  vulgaris  (which 
brought  the  total  volume  to  310  mL)  and,  without  stirring,  the  microalgae  were  allowed 
to  grow  as  lawns.  The  endpoint  for  these  experiments  came  when  the  green  microalgae 
began  to  turn  brown  (beakers:  34  days;  ffasks:  29  days).  Microalgal  yield  was  measured 
as  changes  in  turbidity  (NTU)  using  a  NUL-231  Turbidity  Logger  Sensor  coupled  to  a 
USB-200  USB  Module  (Rochester,  NY,  USA).  No  external  supplies  of  CO2  or  nutrients 
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were  administered  in  these  two  sets  of  cultures  beyond  atmospherically  or  wastewater  sup¬ 
plied  (Table  1 ),  respectively. 

To  measure  and/or  control  the  pH  of  microalgae  cultures,  an  integrated  pH  monitoring 
system  (Doctors  Foster  and  Smith,  Rhinelander,  WI,  USA)  was  used.  The  components 
of  the  system  were  (i)  a  Pintpoint  pH  electrode  (measurement  range:  1.00-14.00;  reso¬ 
lution:  00.00  pH  Unit;  America  Marine  Inc.,  Ridgefield,  CT,  USA),  (ii)  a  pair  of  CO2 
delivery  regulators  (Azoo,  Taipei,  Taiwan)  controlled  by  a  single  solenoid  valve  (Jin  Ben 
Sun  Co.,  Ltd.,  Taoyuan  City,  Taiwan),  and  (iii)  a  Pinpoint  pH  Controller  (America  Ma¬ 
rine  Inc.,  Ridgefield,  CT,  USA)  for  ultimately  powering,  monitoring  and  controlling  all 
the  components  of  the  system.  Before  using  the  automated  pH  monitoring  system  for  ac¬ 
quiring  experimental  values,  a  two-point  pH  electrode  calibration  was  carried  out  using 
a  pH  7  .00  certified  reference  buffer  solution  (to  set  the  lower  detection  limits)  and  a  pH 
10.00  certified  reference  buffer  solution  (to  set  the  slope).  Both  certified  buffer  solutions 
were  supplied  by  Fisher  Chemicals  (Fair  Lawn,  NJ,  USA).  To  monitor  and  regulate  the 
pH  of  30-L  photobioreactors,  the  integrated  system  was  configured  as  shown  in  Fig.  la. 
The  pH  electrode  was  cleaned  and  recalibrated  when  necessary  (usually  once  a  week). 
To  monitor  the  pH  of  the  seven  310-mL  photobioreactor  flasks  (Fig.  lb)  and  beakers 
(not  shown),  the  pH  electrode  was  used  to  gently  stir  the  microalgae  cultures  for  about 
30s  before  taking  measurements.  The  electrode  was  rinsed  with  distilled  water  between 
measurements.  By  design,  the  310-mL  photobioreactors  required  no  additional  CO2  gas 
supplies. 

To  determine  the  empirical  differences  between  an  open  microalgae  photobioreactor 
system  and  a  closed  one  by  turbidity  measurements,  four  separate  30-L  tap  water-based 
C  vulgaris  photobioreactors  similar  to  that  shown  in  Fig.  la  were  set  up  as  follows.  The 
first  two  were  open  to  the  atmosphere  with  one  receiving  additional  supply  of  industrial 
CO2  (Open,  CO2)  and  the  other  receiving  no  additional  CO2  supply  beyond  atmospher¬ 
ically  supplied  (Open,  No  CO2).  The  third  and  the  fourth  photobioreactors  were  closed 
systems  but  similarly  configured  (i.e..  Closed,  CO2  and  Closed,  No  CO2).  Tap  water  (pH 
7.27  zb  0.09,  n  =  4)  was  used  here  because  it  is  the  background  matrix  of  the  wastew¬ 
ater  used  in  this  study.  To  each  photobioreactor  was  added  a  1  L  dilution  of  the  stock 
C.  vulgaris  as  described  earlier  in  this  section,  and  20  mL  of  Guillard’s  (F/2)  Marine  Wa¬ 
ter  Enrichment  Solution  (Carolina  Biological  Supply,  Burlington,  NC,  USA)  every  seven 
days.  Each  photobioreactor  was  equipped  with  a  circulation  pump  (wavemaker.  Fig.  la)  to 
keep  microalgae  in  suspension.  These  four  pumps,  discovered  to  increase  the  temperature 
of  the  microalgae  cultures,  were  turned  on  and  off  automatically  every  30  min  to  keep  the 
temperature  of  the  photobioreactors  to  within  23. 1  ±  0.2°C.  The  light  cycles  and  the  pH  of 
cultures  were  monitored  and  controlled  as  described  earlier  in  this  section.  Turbidity  was 
measured  every  seven  days,  and  on  the  10th  day,  over  21  days. 

Non-linear  modeling  of  turbidity  and  pH  variation  with  time,  and  principal  compo¬ 
nent  analysis  (PCA)  were  performed  using  the  PAST  software  platform  (PAleontolog- 
ical  STatistics;  Hammer  et  al.  2001).  PCA,  an  unsupervised  linear  dimensionality  re¬ 
duction  algorithm  (reviewed  by  Jollife  and  Cadima  2016;  see  also  Gewers  et  al.  2018) 
was  used  to  find  a  more  meaningful  coordinate  system  to  work  with  since  our  vari¬ 
ables  have  different  units  (i.e.,  %WW,  incubation  days  and  pH).  It  was  used  similarly 
to  study  growth  and  biochemical  composition  of  C.  vulgaris  in  different  growth  me¬ 
dia  by  Chia  et  al.  (2013).  Standardizing  data  for  PCA  was  accomplished  by  subtract¬ 
ing  the  mean  value  from  individual  measurement  and  dividing  the  result  by  standard 
deviation. 
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Fig.  2.  The  profiles  of  Clilorella  vulgaris  growth  as  measured  by  the  turbidity  of  the  seven  waste¬ 
water/distilled  water  mixtures  described  in  Fig.  lb.  For  cross  authentication,  duplicate  test  cultures  were 
started  on  different  dates  with  two  different  sets  of  container  geometries  to  run  for  29  days  in  500-mL 
Erlenmeyer  ffasks  (a)  and  34  days  in  500  m-L  beakers  (b).  The  end  of  incubation  was  determined  by  the 
browning  of  the  microalgae  in  the  cultures.  Experimentally  determined  maximum  turbidities  are  shown 
by  the  intersections  of  the  horizontal  dash  lines  and  the  ordinate  in  (c)  and  (d).  Also  shown  at  the  inter¬ 
section  of  each  vertical  dash  line  with  the  abscissa  is  the  optimal  %WW.  The  experimental  coordinate  for 
the  maximum  rabidity  in  the  Erlenmeyer  ffasks  (c)  is  (79%,  333  NTU)  and  in  the  beakers  (d)  is  (78%,  347 
NTU).  Modelling  (not  shown)  suggests  that  the  coordinate  for  maximum  turbidity  should  actually  be  (83%, 
773  NTU). 


Results 

To  (determine  the  optimal  WW/DH2O  ratio  and  incubation  time  required  to  maximize 
biomass  yield,  C.  vulgaris  was  cultured  in  a  pair  of  serially  diluted  wastewater  photobiore¬ 
actors  for  29  and  34  days.  Results  show  that  yield,  as  measured  by  the  turbidity  of  the 
cultures,  is  reproducibly  related  to  time  and  exposed  surface  area  (compare  Erlenmeyer 
flasks  results  displayed  in  Fig.  2a  with  beakers  results  in  Fig.  2b).  The  highest  turbidity 
values  were  measured  after  1 5  days  of  incubation  in  the  83%  WW  culture;  the  lowest  being 
recorded  in  the  100%  WW  irrespective  of  the  length  of  incubation  time. 

To  acquire  potentially  useful  information  for  scaling  up  a  wastewater-based  photobiore¬ 
actor,  a  non-traditional  interpretation  of  our  data  is  provided  in  Figs.  2c-d.  Four  conclu¬ 
sions  can  be  made  from  these  two  plots:  (i)  that  100%  DH2O  or  100%  WW  is  not  the  best 
medium  for  culturing  C.  vulgaris  even  though  100%  DH2O  is  more  accommodating  than 
100%  WW,  (ii)  that  yield  falls  with  increasing  %WW  during  the  first  six  days  of  incuba¬ 
tion,  (iii)  that  80  di  3%  WW  culture  is  the  optimal  range  for  realizing  the  most  biomass  in 
wastewater  provided  incubation  is  allowed  to  continue  for  at  least  29  days  (but  not  more 
than  34  days  because  the  green  microalgae  turns  brown  thereafter  under  this  condition), 
and  (iv)  that  beyond  78-83%  WW,  wastewater  severely  limits  C.  vulgaris  growth. 

To  understand  the  ability  of  C.  vulgaris  to  absorb  atmospheric  CO2,  we  monitored 
changes  in  pH  in  each  of  the  WW/DH2O  cultures  over  time.  The  rationale  was  that  pH  is 
directly  related  to  the  availability  of  dissolved  CO2  in  the  cultures  (Chi  et  al.  201 1)  which 
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Fig.  3.  The  variation  of  pH  with  time  in  the  WW-based  photobioreactors  at  different  %WW:  (a)  500- 
mL  Erlenmeyer  flasks;  (b)  500-mL  beakers.  The  duplicate  experiments  (a)  and  (b)  were  started  at  different 
times  for  cross  authentication,  (c)  Modeling  the  p// variation  in  (a)  and  (b)  to  determine  the  theoretical  pH 
maximum  (i.e.,  p//max  =  ^-69  ±  0.09).  Note  (i)  that  the  p//max  =  ^-69  ±  0.09  is  also  the  mean  of  pH,„ax 
values  in  Table  3,  and  (ii)  that  the  inflection  points  in  (c)  are  near  constant  in  growth  media  containing  0% 
to  33%  WW,  but  vary  widely  between  50%  and  100%  WW. 


can  in  turn  be  directly  linked  to  the  capacity  to  capture  atmospheric  CO2.  Duplicate  results 
from  this  pH  study  are  presented  in  Fig  3a  where  eight  pH  measurements  were  taken  on 
different  days  over  a  period  of  29  days  and  in  Fig.  3b  where  10  measurements  were  taken 
similarly  over  a  period  of  34  days.  In  both  duplicates,  the  lowest  pH  values  were  measured 
in  100%  WW  culture,  and  the  highest  in  the  0%  WW  culture  (DH2O,  pH  7.99  ±  0.06,  n  = 
2).  It  can  also  be  seen  in  both  cases  that  an  upper  pH  ceiling  is  reached  in  each  %WW  cul¬ 
ture  after  10  days  of  incubation.  The  average  pH  value  at  these  ceilings  is  8.7  it  0.1  (/?  =  14). 
This  average  value  reveals  the  enormous  potential  of  a  wastewater-based  photobioreactor 
for  absorbing  CO2. 

Data  from  the  open  WW  photobioreactor  system  (Fig.  1)  show  that  the  relation¬ 
ship  between  turbidity  (x)  and  the  length  of  incubation  in  days  (cl)  in  any  of  our 
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Table  2.  Values  of  constants  D  in  the  Michaelis-Menten  growth  equation  Twhidity  = 

[Unax*(i\/[D  +  (l\  derived  from  regression  analysis  of  turbidity  data  at  each  %WW  (Fig.  2a-b).  They  are 
reported  as  mean±standard  deviation. 


0% 

16% 

33% 

SO'/o 

67% 

83% 

lOO'/o 

323  ±  17 

D  13  ±  0.7 

296  ±  5 

7.9  ±  0.8 

288  ±  18 

10.8  ±  0.6 

379  ±  140 

22.8  ±  15.1 

530  ±  15 
22.5  ±  3.8 

764  ±  13 

41.7  ±  5.5 

239  ±25 

17.6  ±  4.7 

T,nax  is  the  extrapolated  maximum  turbidity. 

x,nax  itnd  D  are  strongly  and  positively  correlated  (Pearson  r  =  0.90). 
Variable  d  is  length  of  incubation  in  days. 


310-mL  wastewater-based  C  vulgaris  photobioreactors  fits  quantitatively  into  the  follow¬ 
ing  (Michaelis-Menten  kinetics)  equation: 


{"^nuix  ^  d) 

{D  -\-  d) 


where  Xnuix^  the  maximum  turbidity,  is  obtained  by  extrapolating  d  in  Eq.  1  beyond  250  days 
in  this  study,  and  D  is  the  number  of  incubation  days  needed  for  x  to  reach  1/2  All  x^ax 

and  D  values  from  this  study  are  listed  in  Table  2. 

To  devise  some  means  of  predicting  pH  changes  during  C.  vulgaris  cultivation,  given  any 
%WW,  a  regression  analysis  (not  shown)  was  performed  for  each  of  the  plots  presented  in 
Fig.  3a-b.  The  results  of  the  regression  analysis  suggest  that  pH  variations  can  be  model 
using  the  logistic  expression  below: 


pf^max 

[1  + 


where  pH„uix  is  the  maximum  pH  to  be  expected  when  d  is  infinitely  large.  The  distribution 
of  pHniax  values  (Table  3)  is  found  to  mirror  that  of  turbidity  (Fig.  2).  This  apparent  corre¬ 
lation  seems  to  suggest  that  pH„uix  is  a  measure  of  the  carbon  requirement  for  C.  vulgaris 
growth  and,  therefore,  an  indirect  measure  of  capacity  to  absorb  CO2  gas.  The  meanings 
of  B  and  C  are  unclear  to  us  at  this  time.  All  values  of  pH,„ax,  B,  and  C  obtained  here  are 
listed  in  Table  3.  Scatter  plots  based  on  this  model  are  presented  in  Fig.  3c  for  comparison 
with  experimental  plots  in  Figs.  3a-b. 

To  determine  the  range  of  %WW  values  and  incubation  days  within  which  photobiore¬ 
actor  setup  such  ours  can  perform  above  background,  PCA  was  applied  to  a  7x  10  correla¬ 
tion  matrix  derived  from  turbidity,  or  pH  values  in  combination  with  corresponding  %WW 
and  number  of  incubation  days.  In  Fig.  4a  is  presented  a  PCA  biplot  for  the  correlation 
matrix  of  turbidity  changes  with  days  of  incubation  in  various  %WW.  The  %WW  markers 


Table  3.  Averages  of  Erlenmeyer  flasks  (Fig.  3a)  and  beakers  (Fig.  3b)  values  for  logistic  constants 
pH  max,  B  and  C  in  modeled  growth  equation  pH  =  pHmax/[^  +  These  values  are  obtained  by  re¬ 

gression  analysis  of  pH  as  a  function  of  incubation  days  at  each  %WW  and  are  reported  as  meanistandard 
deviation. 


O'/o 

16% 

33% 

50'Vi. 

67% 

83% 

100% 

pH  max 

B 

C 

8.82  ±  0.05 
0.09  ±  0.01 
0.15  ±0.06 

8.75  ±  0.00 

0.15  ±0.00 
0.22  ±  0.03 

8.75  ±  0.02 

0.19  ±0.01 
0.20  ±  0.04 

8.68  ±  0.02 

0.23  ±  0.00 
0.19  ±0.00 

8.67  ±  0.09 
0.24  ±0.01 
0.15  ±0.02 

8.61  ±0.19 
0.24  ±  0.04 

0.12  ±  0.04 

8.57  ±  0.08 
0.26  ±0.01 

0.09  ±  0.00 
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Fig.  4.  Scatter  plots  of  two  key  principal  components  (PCI  and  PC2)  from  PCA  to  determine  the  best 
combination  of  %WW,  cultivation  time  (in  days,  d),  and  p//  values  for  optimal  growth,  (a)  The  loading 
of  incubation  days  in  an  ordination  space  illustrating  the  influence  of  %WW  as  vectors  (green  lines)  on 
incubation  days,  (b)  The  scores  of  %WW  in  the  same  space. 


are  shown  as  vectors  (green  lines),  and  the  incubation  days  as  numbers  in‘d’  unit.  There  are 
two  observations  which  can  be  made  from  Fig.  4a.  First,  that  PCI  (which  explains  97.5% 
of  variance)  strongly  correlates  with  seven  of  ten  incubation  days:  0,  3,  6,  24,  27,  32,  and  34 
days.  Of  the  seven  days,  PCI  increases  with  a  set  of  four  incubation  days  longer  than  20  (24, 
27,  32,  and  34  days),  and  decreases  with  a  set  of  three  incubation  days  less  than  10  (0,  3  and 
6  days).  The  former  suggests  that  turbidity  and  the  set  of  days  24,  27,  32,  and  34  (and  days 
in-between)  vary  together.  The  latter  implies  that  incubating  C.  vulgaris  in  wastewater  for 
0,  3  or  6  days  (and  days  most  likely  in-between)  will  lead  to  a  decrease  in  turbidity  as  shown 
here  earlier.  This  axis  (PCI,  Fig.  4a)  can,  therefore,  be  viewed  in  similar  PCA  as  a  measure 
of  the  effect  of  incubation  on  algal  growth;  the  larger  the  value  along  PCI,  the  larger  the 
corresponding  turbidity  value  will  be.  Days  between  10  and  20,  the  ‘Goldilocks’  days,  are 
not  correlated  with  PCI.  Notice  also  that  all  %WW  vectors  (green  lines)  are  accounted 
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Fig.  5.  Determining  the  level  of  gaseous  CO2  required  for  optimal  Clilorella  vulgaris  growth  rate  in 
distilled  H2O  (30  L,  Fig.  la).  For  the  first  21  days  (left  of  vertical  dashed  line),  Clilorella  vulgaris  was  supplied 
with  nutrients  as  described  in  the  text,  but  deprived  of  nutrients  thereafter.  Temperature  was  maintained  at 
23.1  ±  0.2°C  over  the  42-days  of  experiment. 


for  along  PCI  with  positive  loadings.  This  simply  restates  the  benefit  derives  from  having 
wastewater  in  C  vulgaris  cultures  since  the  larger  the  values  are  along  PCI,  the  better. 

PCI,  however,  fails  to  offer  the  %WW  range  that  would  support  optimal  C.  vulgaris 
growth.  This  information  is  provided  in  the  loadings  of  WW  vectors  along  the  second 
principal  component,  PC2.  It  can  be  seen  that  PC2,  which  explains  1.8%  of  variation, 
correlates  positively  with  the  0,  16,  33,  and  50%  WW  vectors,  but  negatively,  with  the  67, 
83,  and  100%  WW  vectors.  The  0,  16,  33,  and  50%  WW  cultures  were  shown  earlier  in 
this  section  to  lead  to  less  than  optimal  turbidity  values  and  hence  do  not  support  high 
C.  vulgaris  yield.  The  83%  WW  culture  on  the  other  hand  was  shown  to  be  the  optimal 
WW  level  suitable  for  growth  (followed  by  the  67%  WW  culture;  Fig.  2).  Based  on  these 
observations,  values  along  PC2  can  therefore  be  considered  a  measure  of  the  capacity  of 
wastewater  to  support  C.  vulgaris  growth;  the  small  the  value,  the  large  the  turbidity  will 
be.  The  83%  WW  then  stands  out  as  the  best  growth  medium  for  C.  vulgaris  since  it  has 
the  lowest  value  along  PC2. 

Fig.  4b  is  a  PC  A  scatter  plot  for  pH.  It  can  be  seen  that  PCI  (which  explains  92.3% 
of  variance)  is  strongly  and  positively  correlated  with  lower  %WW  (0  and  16%  WW),  but 
negatively  correlated  with  higher  %WW  values  (83  and  100%).  Because  evidence  shows 
that  lower  %WW  does  not  support  growth  optimally,  but  83%  is,  this  principal  component 
can  be  used  as  a  measure  of  the  best  condition  for  absorbing  CO2  (since  the  p//  of  an  open 
photobioreactor  is  closely  linked  to  the  availability  of  CO2  in  the  atmosphere).  The  smaller 
the  values  along  PCI,  the  higher  the  capacity  of  the  microalgae  to  absorb  CO2  gas.  Note 
that  the  clustering  of  wastewater  fractions  in  Fig.  4b  mirrors  the  distribution  of  B  values 
with  %WW  in  Table  3  which  may  be  suggesting  the  meaning  of  constant  B  in  Eq.  2. 

To  understand  the  role  of  an  open  or  a  closed  system  against  as  far  as  CO2  supply  to  a 
C.  vulgaris  culture  is  concerned,  four  separate  30-L  tap  water-based  photobioreactors  were 
set  up.  Fig.  5  is  a  summary  of  the  results  of  these  experiments.  It  can  be  seen  that  growth 
was  generally  exponential  in  each  photobioreactor  over  the  21  days  of  experiment  (left  side 
of  dashed  line  in  Fig.  5).  During  this  time,  the  open  system  registered  the  best  performance 
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with  additional  CO2.  The  worse  performance  was  observed  in  the  two  photobioreactors 
configured  to  be  open  with  no  additional  CO2  supply  or  closed  with  additional  supply  of 
CO2.  To  reaffirm  the  importance  of  non-maternal  nutrients  to  C  vulgaris,  the  four  photo¬ 
bioreactors  above  were  incubated  for  19  more  days  without  added  nutrients.  The  results, 
summarized  on  the  right  side  of  Fig.  5  (days  between  21  and  40)  show  that  biomass  pro¬ 
duction  declined  significantly  during  the  period,  even  in  the  best  photobioreactor  configu¬ 
ration  (i.e.,  open  with  additional  CO2  supply). 


Discussion 

We  have  shown  that  the  temporal  growth  profile  of  C.  vulgaris  in  a  wastewater-based 
photobioreactor  follows  a  nonlinear  kinetics  (Figs.  2a-b).  This  is  consistent  with  what 
others  have  observed  (Ammar  2016).  However,  that  C  vulgaris  exhibits  clearly  seg¬ 
mented  growth  response  to  wastewater  as  shown  in  Figs.  2c-d  is  a  novelty.  In  this  newly 
discovered  profile,  the  accumulation  of  C.  vulgaris  biomass  declines  progressively  in  low 
levels  of  wastewater  (0-50%  WW),  recovers  rapid  in  the  50-83%  WW  range  and  then  falls 
precipitously  in  the  83%  to  100%  WW  range.  To  explain  this  profile,  it  is  useful  to  think  of 
the  decline  in  low  levels  of  wastewater  as  a  reflection  of  C.  vulgaris  adapting  to  wastewater- 
induced  shock  (Wang  et  al.  2009).  This  is  justifiable  because  the  burden  of  adjustment 
logically  increases  with  the  proportion  of  wastewater  in  cultures.  The  rapid  biomass  in¬ 
crease  in  the  50-83%  WW  range  then  can  be  thought  of  as  the  benefits  microalgae  derive 
from  wastewater,  however  miniscule  they  may  be  in  low  %WW  cultures  is  enough  to  over¬ 
come  the  wastewater-related  shock  in  this  %WW  range.  Beyond  this  %WW  range,  the 
precipitous  drop  in  biomass  production  may  be  explained  hypothetically  as  a  breach  of 
microalgae  cellular  biochemistry  by  high  levels  of  wastewater  components.  More  studies 
are  required  to  test  these  hypotheses. 

The  role  that  ‘maternal’  nutrients  play  in  growth,  which  is  clearly  exhibited  in  Fig.  2 
(0%  WW  cultures),  is  worth  mentioning  here.  The  fact  that  C.  vulgaris  grows  at  all  and 
more  so  in  100%  DH2O  means  that  maternally  stored  nutrients  (in  the  cytoplasm  and 
in  the  inter-thylakoid  regions  of  C.  vulgaris  chloroplast;  Safi  et  al.  2014)  are  used,  we  be¬ 
lieve,  for  autonomous  growth.  This  continues  until  the  maternal  nutrients  are  exhausted  (in 
about  the  sixth  day  of  incubation;  Fig.  2d)  whereupon  wastewater  becomes  the  source  of 
nutrients. 

This  study  finds  that  the  optimal  conditions  for  culturing  C  vulgaris  in  an  open 
wastewater-based  photobioreactor  is  80  ±  3%  WW  over  a  29-34  incubation  days  window. 
It  also  finds  that  additional  metabolic  carbon  source  beyond  atmospheric  CO2  is  required 
to  maximize  the  biomass  production  in  a  photobioreactor  (i.e..  Open,  CO2  system.  Fig.  5). 
The  importance  of  an  open  photobioreactor  system  also  needs  to  be  emphasized  here  for 
two  reasons.  First,  studies  show  that,  without  additional  source  of  carbon  (industrial  CO2 
for  example),  the  capacity  of  C.  vulgaris  cells  to  photosynthesize  falls  which  in  turn  causes 
a  drop  in  biomass  production  (Lundquist  et  al.  2009;  Gebreslassie  2013;  Raven  2017).  In¬ 
deed,  most  photobioreactors  fail  because  of  inadequate  supplies  of  metabolic  carbon  and 
yet  one  cheap  source  of  CO2,  the  flue  gas  from  power  plants,  is  always  available  (Gebres¬ 
lassie,  2013;  Kumar  et  al.  2011).  Secondly,  O2  gas,  the  byproduct  of  photosynthesis,  is  toxic 
to  algae  (Molina-Grima  et  al.  2001).  It  has  been  shown  that  microalgae  can  only  tolerate 
up  to  400%  of  O2  above  the  level  in  O2  saturated  air  (Chisti  2007).  Having  an  open  wastew¬ 
ater  photobioreactor  would  therefore  reduce  toxicity  by  allowing  O2  produced  during 
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photosynthesis  to  escape  (while  still  leaving  enough  for  microbes  to  use  in  the  breaking 
down  of  organic  matter  in  wastewater;  Azov  et  al.  1982;  Abdel-Raouf  et  al.  2012). 

This  study  identified  five  elements  of  a  wastewater-based  C  vulgaris  photobioreactor 
pHnutx-:  D,  B,  C)  which  could  be  useful  for  forecasting  productivity  even  though  the 
meanings  of  B  and  C  is  still  to  be  determined.  For  any  chosen  %WW,  we  have  introduced 
(i)  the  concept  of  the  minimum  days  (D)  a  static  photobioreactor  would  need  to  produce 
biomass  above  background,  (ii)  the  concept  of  pH„,ax,  the  upper  ceiling  of  pH  which  can  be 
used  to  predict  the  capacity  of  a  wastewater-based  C.  vulgaris  photobioreactor  to  absorb 
CO2  in  an  open  photobioreactor  system,  and  (iii)  the  concept  of  T,„ax^  the  upper  ceiling  of 
biomass  production  when  incubation  is  continued  indefinitely  beyond  the  minimum  days, 
D  (i.e.,  cl  >>  D;  Eq.  1).  One  property  of  needs  elaboration.  Unlike  pH,„ax  and  D, 
the  values  of  x„iax  cannot  be  obtained  directly  by  C.  vulgaris  cultivation  because  of  the 
limits  imposed  on  production  by  the  carrying  capacity  of  photobioreactors  and  the  long 
incubation  time  required  to  achieve  it.  For  example,  the  average  carrying  capacity  of  our 
310-mL  system  was  340  NTU  (the  mean  of  333  NTU  and  347  NTU;  Fig.  2a  and  Fig.  2b, 
respectively)  and  yet  was  determined  to  be  773  NTU  in  a  83%  WW  culture  (after 
550  years  of  incubation). 

A  number  of  factors  were  confirmed  in  this  study  to  have  significant  consequence  on  the 
effectiveness  of  using  microalgae  to  polish  municipal  wastewater.  For  example,  we  find  that 
efficient  removal  of  nutrients  from  wastewater  as  measured  by  increase  in  turbidity  depends 
on  constant  agitation  of  the  C.  vulgaris  cultures.  This  is  of  consequence  in  municipalities 
where  microalgae  have  been  adopted  to  polish  wastewater,  but  agitation  has  not  been  incor¬ 
porated  into  their  photobioreactor  designs  (Craggs  et  al.  2003).  Exacerbating  this  problem 
is  the  fact  that  the  C:N  ratio  in  a  typical  municipal  wastewater  is  low  when  compared  to 
the  carbon  requirements  of  microalgal  biomass.  This  low  C:N  ratio  limits  the  capacity  of 
microalgae  to  produce  biomass  in  wastewater  (Benemann  2003).  This  imbalance  is  usually 
manifested  as  elevated  pH  levels  during  incubation  (resulting  from  the  use  of  bicarbonate 
ions  as  a  CO2  source  for  algal  photosynthesis,  releasing  hydroxide  ions;  Craggs  et  al.  201 1) 
which  in  turn  prevents  growth  of  both  microalgae  and  bacteria  which  degrade  the  organic 
compounds  in  wastewater  (Azov  et  al.  1982).  This  scenario  was  observed  in  photobioreac¬ 
tors  which  were  not  supplied  with  additional  CO2  (Figs.  3,  5).  Boosting  carbon  supply  to  a 
wastewater-based  culture  would  not  only  enhance  microalgal  yield  as  has  been  shown  here, 
it  would  also  fulfill  a  goal  to  scrub  CO2  gas  from  sources  beyond  the  atmosphere.  Note  that 
the  main  function  of  microalgae  for  polishing  municipal  wastewater  is  to  remove  N,  and 
not  P  (Choi  and  Lee  2015).  This  is  because  the  large  N:P  ratio  requirement  for  biomass 
production  in  microalgae  means  no  additional  microalgal  biomass  production  is  needed 
above  that  required  to  assimilate  N  to  remove  P  from  municipal  wastewater  (Lau  et  al. 
1995;  Benemann  2003;  Craggs  et  al.  2011). 

The  environmental  remediation  approach  proposed  in  this  study  will  results  in  mas¬ 
sive  amounts  of  algal  biomass,  the  fate  of  which  must  be  addressed  appropriately.  Fortu¬ 
nately,  some  means  for  doing  so  already  exist  even  though  they  were  not  designed  with  ad¬ 
dressing  environmental  issues  in  mind.  For  example,  microalgae  biomass  has  been  turned 
into  chemical  feedstock  for  biodegradable  plastics  (Lambert  and  Wagner  2017),  lubricants 
(Ruiz  et  al  2016),  fertilizers  (Uysal  et  al.  2015),  and  pharmaceutical/nutraceutical  products 
(Jha  et  al.  2017).  In  the  skin-care  natural  product  industry,  the  capacity  of  microalgae  to 
developed  sun  blocking  agents  to  protect  and  heal  themselves  from  the  damaging  impacts 
of  exposure  to  solar  radiation  and  other  environmental  hazards  has  been  harnessed  to  pro¬ 
tect  human  skin  (Stolz  and  Obermayer  2005).  In  some  parts  of  the  world,  microalgae  have 
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been  used  sustainably  as  sources  of  animal  feed  because  of  their  ability  to  concentrate  car¬ 
bohydrates,  proteins  and  vegetable  oils,  micronutrients,  vitamins,  and  valuable  pigments 
(Lambert  and  Wagner  2017).  In  an  application  incompatible  with  our  goal,  microalgae 
biomass  has  also  been  used  as  human  dietary  supplements  (Yaakob  et  al.  2014).  And  of  re¬ 
cent,  a  technology  has  emerged  that  has  potential  to  convert  microalgal  biomass  into  a  host 
of  products.  Currently,  this  proprietary  technology  employs  high  energy  electron  beams  to 
convert  cellulose  in  corn  cobs  to  edible  sugar,  alcohol,  biodegradable  plastics  and  several 
other  products  on  industrial  scaled  Microalgal  biomass  could  be  treated  similarly. 

Capturing  and  storing  CO2  emission  is  one  feasibly  way  of  mitigating  climate  change. 
How  the  captured  CO2  is  stored  however  will  have  an  impact  on  the  effectiveness  of 
long-term  sequestration.  For  example,  the  ocean’s  capacity  to  absorb  and  store  huge  quan¬ 
tities  of  CO2  could  be  used  to  sequester  CO2.  This  though  comes  with  environmental  cost 
and  legal  problems.  First,  several  studies  have  shown  that  sequestering  CO2  in  the  ocean 
leads  to  global  warming  which  ultimately  leads  to  excessive  plant  growth.  As  pointed  out  in 
the  Introduction  section,  this  excessive  plant  growth  creates  dead  zones  in  the  ocean  (Shaf¬ 
fer  2010).  Secondly,  the  Clean  Water  Act  makes  it  impossible  to  sequester  CO2  into  the 
ocean. The  Act  prohibits  point  source  discharges  into  navigable  waters  in  the  USA  with¬ 
out  a  permit;  it  also  ensures  that  marine  environments  will  not  suffer  unreasonable  degra¬ 
dation  or  irreparable  harm  from  anthropogenic  activities.  One  alternative  to  the  ocean  is 
geological  sequestration  (Shaffer  2010).  This  alternative  involves  injecting  captured  CO2 
into  geologic  environments  (Aarnes  et  al.,  2010;  Plasynski  et  al.,  2011).  However,  the  alter¬ 
native  also  comes  with  one  major  environmental  concern,  that  of  leaks  from  storage  sites. 


Conclusions 

We  have  demonstrated  that  a  microalgae-based  photobioreactor  can  offer  a  potential  so¬ 
lution  to  lowering  the  level  of  atmospheric  CO2  while  cleaning  nutrients-laden  municipal 
wastewater  at  the  same  time.  By  culturing  C.  vulgaris  in  various  photobioreactor  systems 
containing  wastewater,  we  have  shown  that  wastewater  plays  two  competing  roles  during 
C  vulgaris  growth:  a  mostly  detrimental  and  dominant  role  at  low  levels  where  we  believe 
wastewater  interferes  with  the  metabolic  mobilization  of  the  maternal  nutrients  required 
for  the  initial  growth,  and  a  mostly  supportive  role  that  is  very  evident  in  the  50-83%  WW 
range  cultures.  As  such,  the  optimal  conditions  for  C  vulgaris  to  absorb  atmospheric  CO2 
are  found  to  be  78-83%  dilution  of  wastewater  (the  baseline  is  33%  WW)  and  24  days 
of  incubation  in  open  air.  We  have  also  shown  that  principal  components  analysis  can 
be  applied  to  rudimentarily  acquired  data  such  as  ours  to  aid  in  determining  the  opti¬ 
mal  wastewater  requirement  of  a  photobioreactor.  The  size  of  the  chosen  photobioreactor 
would  be  limited  by  the  carrying  capacity  of  photobioreactor. 
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Abstract. — The  eastern  gray  squirrel,  Sciurus  carolinensis  (EGS)  has  been  introduced 
to  California  and  has  expanded  its  geographic  range  since  initial  introductions.  In  this 
study  we  projected  the  potential  future  geographic  range  of  the  EGS  in  California  us¬ 
ing  Maxent  to  create  an  ecological  niche  model.  Location  data  were  obtained  over  the 
time  period  of  2004-2015  from  museum  specimens,  wildlife  rehabilitation  centers,  the 
California  Department  of  Public  Health,  the  California  Roadkill  Observation  System, 
and  non-iNaturalist  citizen  science  observations.  Research  grade  data  from  iNatural- 
ist  was  obtained  over  the  time  period  of  2004-2018.  Range  and  habitat  suitability 
maps  were  developed  by  mapping  in  ArcGIS.  Three  threshold  selection  methods  were 
used  to  create  different  estimates  of  the  potential  future  range  of  the  EGS  in  Cali¬ 
fornia.  The  first  method  used  the  10th  percentile  logistic  threshold,  the  second  used 
the  minimum  training  presence  logistic  threshold,  and  the  third  used  Jenks  Natural 
Breaks.  We  propose  that  Jenks  Natural  Breaks  has  distinct  advantages  over  the  other 
two  methods  for  estimating  the  potential  future  range  of  the  introduced  EGS  in  Cal¬ 
ifornia,  because  it  provides  information  on  the  habitat  suitability  ranking  throughout 
California,  whereas  the  other  methods  only  provide  a  binary  suitable/ unsuitable  map. 


The  objective  of  this  study  was  to  develop  an  ecological  niche  model  (ENM)  with  an  ap¬ 
propriate  threshold  value  that  could  best  identify  potential  range  expansion  of  the  invasive 
eastern  gray  squirrel,  Sciurus  carolinensis  (EGS)  within  California  and  in  the  future,  project 
potential  areas  of  overlap  with  congeners.  The  EGS  is  native  to  the  deciduous  forests  of 
the  eastern  United  States  (Koprowski  1994).  The  species  was  introduced  to  California  in 
1939  at  Stanford  University  and  in  1943  at  Golden  Gate  Park  in  San  Francisco,  but  it 
may  have  been  introduced  earlier  by  settlers  who  frequently  introduced  species  from  their 
homes  in  the  eastern  United  States  (Byrne  1979).  Populations  of  the  EGS  currently  exist  in 
developed  and  forested  areas  of  California.  They  are  currently  widespread  throughout  cen¬ 
tral  California,  with  concentrations  around  Sacramento,  both  peninsulas  of  San  Erancisco 
Bay,  and  San  Jose,  and  smaller  populations  in  the  Central  Valley.  Additionally,  populations 
are  spreading  from  Santa  Cruz  into  the  Santa  Cruz  Mountains  and  Monterey  Peninsula 
(Creley  and  Muchlinski  2017). 

Sciurus  carolinensis  carolinensis  is  most  likely  the  major  subspecies  present  in  Califor¬ 
nia,  as  determined  by  the  coat  color  and  physical  characteristics  of  observed  squirrels. 


*  Corresponding  Author:  carlycreley@gmail.com 
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S.  c.  cawlmensis  has  a  gray  dorsum  with  a  cinnamon  wash  sometimes  present  on  the 
dorsum  and  hips.  The  tail  is  the  same  shade  of  gray,  with  a  light  white  frosting  on  the 
tips  of  the  hairs.  A  white  eye  ring  is  usually  visible  (Thorington  et  al.  2012).  Some  EGSs 
in  California  are  melanistic,  a  common  trait  in  the  northern  portion  of  the  native  range 
(Thorington  et  al.  2012).  The  EGS  has  a  broad  diet  (Bertolino  2008),  can  establish  a  pop¬ 
ulation  from  a  small  number  of  founders  (Wood  et  al.  2007),  can  survive  and  reproduce 
in  urban,  suburban,  or  natural  habitats,  and  has  a  favorable  public  perception  (Bertolino 
and  Genovesi  2003),  which  all  contribute  to  its  invasive  success.  The  species  has  been  intro¬ 
duced  to  the  western  United  States,  Europe,  Africa,  and  Australia  (Bertolino  2008;  Benson 
2013;  Bertolino  and  Lurz  2013;  Bertolino  et  al.  2014).  Populations  have  been  associated 
with  negative  ecosystem  effects,  the  decline  of  native  species,  and  damage  to  forests  in  the 
United  Kingdom,  Ireland,  Italy,  and  parts  of  western  North  America  (Gurnell  et  al.  2004; 
Bertolino  2008;  Benson  2013;  Bertolino  and  Lurz  2013). 

The  original  native  range  of  the  EGS  consists  of  mature,  continuous  woodlands  over 
40  ha  in  size,  with  diverse  woody  understories  and  tree  species  such  as  oak  {Quercus), 
hickory  {Ccirya),  and  walnut  {Juglans)  (Koprowski  1994).  However,  EGSs  can  also  live 
in  urban  and  suburban  environments,  even  with  relatively  few  mature  trees  (Thorington 
et  al.  2012).  EGSs  move  primarily  along  river  corridors,  secondarily  on  roads/right-of- 
ways,  and  thirdly  on  tracks/paths  (Stevenson  et  al.,  2013).  Their  ability  to  live  in  developed 
environments  may  significantly  increase  their  dispersal  capability. 

In  order  to  understand  the  future  potential  distribution  of  the  EGS,  it  is  critical  to  model 
their  potential  landscape  and  ecological  niche  occupancy.  Creley  and  Muchlinski  (2017) 
mapped  the  species  distribution  within  California  as  of  2015,  but  ENMs  had  not  been 
made.  The  methodologies  for  creating  ENMs  for  invasive  species  are  not  well  established 
(Aguierre-Gutierrez  et  al.  2013;  Uden  et  al.  2015).  It  was  important  to  map  a  range  of  suit¬ 
able  habitat  estimates  in  order  to  prevent  drawing  conclusions  from  one  arbitrary  thresh¬ 
old,  as  cautioned  against  by  Merow  et  al.  (2013).  Standard  thresholds  that  have  been  used 
in  previous  studies  to  produce  ENMs  include  the  10^^  percentile  logistic  threshold  (Belar- 
main  Fandohan  et  al.  2015;  Chalghaf  et  al.  2016;  Sage  et  al.  2017),  which  could  produce 
a  very  conservative  estimate  of  potential  future  range  for  an  invasive  species,  and  the  min¬ 
imum  training  presence  logistic  threshold  (Beane  et  al.  2013;  Coudrat  and  Nekaris  2013; 
Botero-Delgadillo  et  al.  2015),  which  could  produce  an  overestimate  of  potential  future 
range.  We  used  these  two  standard  thresholds  in  the  present  study,  as  well  as  a  third  method 
-  selecting  a  threshold  value  based  on  Jenks  Natural  Breaks  (Jenks  1967).  Previously,  Col- 
nar  and  Landis  (2007)  developed  a  regional  risk  assessment  for  the  European  green  crab, 
Carcinus  inaenas,  at  Cherry  Point,  Washington,  USA  using  Jenks  Natural  Breaks,  and 
Schleier  III  and  Sing  (2008)  used  it  to  partition  an  overall  risk  score  for  the  introduction  of 
Gahusia  ajfinis  (western  mosquitofish)  into  Montana  watersheds.  Beyond  invasive  species 
modeling,  Jenks  Natural  Breaks  have  been  used  to  classify  groundwater  into  zones  of  vul¬ 
nerability  for  nitrogen  contamination  in  Florida’s  aquifers  (Cui  et  al.  2016),  to  rank  the 
susceptibility  of  locations  to  terrorist  actions  (Patterson  and  Apostolakis  2007),  and  to  as¬ 
sess  the  risk  of  flooding  in  the  Bengawan  Solo  River  basin  in  Indonesia  (Rahadianto  et  al. 
2015). 

The  model  used  location  and  environmental  data  from  the  invaded  range  in  California 
and  the  native  range  for  S.  cawlmensis  in  the  eastern  United  States,  which  encompasses 
a  wider  range  of  environmental  conditions  than  those  found  within  the  current  range  in 
California,  to  estimate  habitat  suitability.  Our  methods  may  be  applicable  to  further  stud¬ 
ies  on  invasive  species  modeling,  for  which  methods  of  estimating  the  potential  range,  as 
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opposed  to  the  current  range,  are  limited  and  have  not  been  consistently  evaluated 
(Aguierre-Gutierrez  et  al.  2013;  Uden  et  al.  2015). 


Materials  and  Methods 

We  obtained  presence  only  location  data  for  2004  to  2018.  The  model  included  3,725 
spatial  location  points  of  the  EGS  in  California,  as  well  as  8,988  points  from  across  the 
native  range  in  the  eastern  United  States.  We  obtained  presence  points  in  California  from 
iDigBio,  the  Global  Biodiversity  Information  Facility  (GBIF),  Vertnet,  wildlife  rehabili¬ 
tation  centers,  the  California  Department  of  Public  Health’s  West  Nile  Virus  Surveillance 
Program,  the  California  Roadkill  Observation  System  operated  by  the  University  of  Cali¬ 
fornia,  Davis  (Waetjen  and  Shilling  2017),  iNaturalist,  and  the  authors.  Focation  data  from 
the  native  range  of  the  FGS  across  the  United  States  are  from  iNaturalist.  We  filtered  the 
data  to  include  only  the  native  range  of  S.  carolinensis,  according  to  Koprowski  (1994).  We 
excluded  regions  from  which  the  FGS  is  nonnative,  and  areas  outside  of  the  United  States. 

VertNet,  a  National  Science  Foundation  funded  project,  makes  museum-curated 
biodiversity-data  free  and  available  on  the  web,  while  the  Global  Biodiversity  Informa¬ 
tion  Facility  provides  open  international  data.  The  iNaturalist  sightings  were  filtered  to 
include  only  those  that  were  open  access  and  research  grade,  which  include  an  observation 
date,  photo,  coordinates,  and  in  which  the  species  identification  has  been  verified  by  at  least 
one  other  user.  Records  in  biodiversity  databases  are  constantly  changed  and  updated,  so 
all  data  from  iDigBio,  GBIF,  and  Vertnet  were  obtained  on  31  August  2015.  Records  from 
iNaturalist  were  obtained  on  12  October  2018.  Some  redundancy  may  exist  between  the 
databases,  but  Maxent  automatically  removes  replicates,  so  these  duplicate  observations 
did  not  change  the  projections. 

Since  all  sources  are  likely  to  include  some  misidentifications  of  related  species  identified 
as  FGS,  reports  from  outside  of  the  previously  published  range  were  scrutinized  for  accu¬ 
racy.  Field  surveys  were  conducted  in  regions  that  had  not  been  included  in  prior  range 
maps  but  that  had  numerous  reports,  including  the  Santa  Cruz  Mountains,  the  Central 
Valley,  and  southern  California.  Records  that  could  not  be  corroborated  were  expunged. 
In  the  case  of  the  California  Roadkill  Observation  System  data,  82%  of  FGS  with  images 
were  correctly  identified.  The  remainders  were  misidentified  as  California  ground  squirrel, 
Douglas  squirrel,  or  western  gray  squirrel  (WGS). 

Data  used  in  this  manuscript  are  available  for  use  by  others  under  a  Creative  Com¬ 
mons  By  Attribution  Non-Commercial  4.0  International  Ficense.  Observations  ob¬ 
tained  through  iNaturalist  are  utilized  under  a  Creative  Commons  By  Attribution  Non- 
Commercial  Ficense  or  from  observations  that  are  in  the  Public  Domain.  The  names  of 
GBIF  and  iNaturalist  contributors  can  be  obtained  through  a  search  of  the  California 
location  data  using  the  data  posted  at  doi:  10.1 3 140/RG.2. 2. 24275. 84004  and  the  native 
range  of  the  United  States  data  at  doi:  10. 13140/RG.2. 2. 17564.95360. 

We  converted  location  data  into  geographic  coordinates  with  Google  Maps.  We  spa¬ 
tially  rarefied  the  data  to  one  point  per  51.8  ha  (0.25  mi“)  using  the  Spatially  Rarefy  tool 
in  the  SDM  toolbox  (Brown  2014)  in  order  to  eliminate  spatial  clusters  that  could  cause 
the  model  to  be  overfit  to  the  environmental  biases  of  those  points  (Boria  et  al.  2014).  The 
rarefied  presence  data  were  reduced  from  12,713  to  5,627  points.  Since  the  predictive  accu¬ 
racy  of  the  maps  and  the  ability  of  the  models  to  project  presence  data  was  more  important 
than  identifying  the  tolerance  ranges  of  the  species,  we  did  not  remove  highly  correlated 
variables  (Merow  et  al.  2013). 


RANGE  EXPANSION  OF  THE  EASTERN  GRAY  SQUIRREL 


61 


For  the  environmental  background  we  used  the  bio  1  through  bio  19  variables  from  BIO- 
CLIMl  remote  sensing  data  (Hijmans  et  al.  2005),  which  include  quarterly  and  annual 
temperature  and  precipitation  trends,  with  other  annual  environmental  trends.  Addition¬ 
ally,  we  used  monthly  precipitation,  monthly  maximum  temperature,  monthly  minimum 
temperature,  altitude,  impervious  surface,  land  cover,  and  tree  canopy  from  the  United 
States  Geological  Survey  (USGS  2016).  The  full  list  of  biotic  and  abiotic  factors  used  in 
our  model  is  available  at  doi:  10.1 3 140/RG.2.2. 15802.70085.  We  paired  the  location  pres¬ 
ence  data  with  environmental  background  data  throughout  the  contiguous  United  States. 
The  model  can  predict  habitat  suitability  for  areas  that  the  species  is  equally  likely  to  reach 
(Merow  et  al.  2013).  The  environmental  background  data  covers  all  reasonable  possible 
distributions  (Saupe  et  al.  2012). 

We  selected  Maxent  to  create  the  ENM  because  of  its  high  performance  at  estimating 
local  occurrences  with  small  to  medium  sample  sizes  (Elith  et  al.  2006;  Aguierre-Gutierrez 
et  al.  2013;  Ng  and  Jorda,  2001).  In  order  to  allow  the  model  to  reach  the  default  0.00001 
convergence  level,  we  allowed  a  maximum  of  5,000  iterations  (Young  et  al.  2011),  set  the 
regularization  parameter  to  the  default  value  of  one  in  order  to  reduce  overfitting  (Merow 
et  al.  2013),  and  created  15  replicates  for  the  model  using  subsampling  (Young  et  al.  2011). 
We  set  aside  25%  of  the  data  for  testing,  and  used  75%  for  training  the  model.  We  used  the 
random  seed  option  in  order  to  increase  the  randomness  of  the  runs  (Jobe  and  Zank  2006; 
Young  et  al.  2011).  We  adjusted  the  sample  radius  to  -100,  and  did  not  extrapolate.  The 
logistic  output  using  the  default  x  value  of  0.5  was  selected  because  the  actual  probability 
of  the  EGS  being  present  in  suitable  habitat  is  unknown. 

We  produced  habitat  suitability  maps  by  importing  the  rasters  for  the  average  of  the 
fifteen  replicates  for  each  model  into  ESRTs  ArcMap  10.3.1  (ArcGIS®  10.3.1;  Esri  soft¬ 
ware)  using  the  NAD  1983  California  (Teale)  Albers  (Meters)  projected  coordinate  system 
and  the  GCS  North  American  1983  (NAD  1983)  datum.  When  coordinates  were  provided 
without  a  datum,  they  were  assumed  to  be  in  NAD83.  Each  pixel  encompasses  approxi¬ 
mately  720  square  meters.  We  clipped  the  output  rasters  to  the  shape  of  California,  using 
the  United  States  Census  Bureau’s  Tiger/Line  2010  (United  States  Census  Bureau  2010). 
We  used  topographic  and  political  basemap  layers  from  ESRI. 

We  used  three  thresholds  to  create  potential  range  estimates.  The  first  estimate  of  the  po¬ 
tential  range  was  established  by  using  the  10th  percentile  training  presence  logistic  thresh¬ 
old,  which  set  the  threshold  at  the  level  where  90  percent  of  the  presence  points  were  in 
raster  squares  with  at  least  the  threshold  score.  We  created  a  second  estimate  of  the  po¬ 
tential  range  using  the  minimum  training  presence  logistic  threshold,  which  is  set  at  the 
level  of  the  lowest  scoring  occupied  square.  We  manually  classified  the  resulting  raster  of 
California  for  the  first  two  thresholds  by  setting  the  upper  bound  of  the  unsuitable  habi¬ 
tat  class  at  the  threshold  level  for  each  map,  as  indicated  in  the  Results  table  created  by 
Maxent. 

Third,  we  used  Jenks  Natural  Breaks  to  group  the  resulting  raster  into  five  classes 
based  on  natural  breaks  in  the  data.  The  approach  grouped  the  relative  habitat  suitabil¬ 
ity  ranking  of  each  raster  square  (Rahadianto  et  al.  2015)  by  similar  values  and  max¬ 
imized  the  difference  between  the  classes  (ESRI  2016).  We  created  5  habitat  categories 
(HC)  from  the  groupings  of  data.  HCl  represents  where  the  species  is  now  found  and 
adjacent  highly  suitable  habitat,  HC2,  HC3,  and  HC4  represent  decreasingly  suitable, 
but  still  suitable  habitat,  while  HC5  represents  unsuitable  habitat.  Therefore,  the  thresh¬ 
old  for  suitable  habitat  was  set  at  the  lowest  relative  habitat  suitability  ranking  score 
of  HC4. 
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Results 

The  high  value  of  the  Test  Area  Under  the  Curve  ( AUC)  of  the  Receiver  Operating  Char¬ 
acteristic  (ROC)  Plot  (0.8272  on  a  scale  of  0.5  =  random  association  of  data  to  model,  to 
1.0  =  perfect  association  of  data  to  the  model)  and  Training  AUC  (0.828)  indicate  an 
excellent  fit  of  the  data  to  the  model  produced  by  Maxent.  The  plot  is  available  at  doi: 
1 0. 1 3 140/RG.2. 2. 1253 1.78883.  However,  since  the  ROC  Plot  is  constructed  from  all  pos¬ 
sible  threshold  values  and  the  plot  does  not  yield  information  about  any  specific  threshold 
value,  the  selection  of  a  biologically  meaningful  threshold  value  for  graphing  results  of  the 
model  is  critical  for  obtaining  the  most  accurate  map  of  potential  future  range. 

Three  dififerent  threshold  values  were  used  for  the  range  of  suitable  to  unsuitable  habitat 
found  using  the  ENM  for  the  EGS  in  California.  Use  of  the  10^*^  percentile  logistic  thresh¬ 
old  (Fig.  1)  produced  the  most  limited  estimate  of  potential  future  range.  In  the  binary 
map,  suitable  habitat  scored  above  the  threshold  value  of  0.3051  (Table  1).  The  potential 
range  produced  using  this  threshold  value  is  very  similar  to  the  current  range  of  the  species 
in  northern  and  central  California  as  of  2018.  Therefore,  the  time  frame  of  projected  range 
expansion  is  very  limited.  EGSs  currently  exist  outside  of  habitat  that  is  predicted  suitable 
through  use  of  this  threshold,  as  the  method  excludes  the  10  percent  of  observation  lo¬ 
cations  with  the  lowest  relative  scores  in  order  to  produce  the  threshold  value.  EGSs  are 
currently  found  outside  of  projected  suitable  habitat  near  Salinas,  Tracy,  Modesto,  and  east 
and  south  of  Sacramento  along  the  Sierra  Nevada  foothills  from  north  of  the  American 
River  to  Columbia. 

Use  of  the  minimum  training  presence  logistic  threshold  (Fig.  2)  also  produced  a  binary 
output  map,  but  in  this  case  the  logistic  threshold  value  was  reduced  to  0.0299  (Table  1).  A 
very  large  increase  in  projected  suitable  habitat  for  the  EGS  in  California  results  from  the 
reduction  in  threshold  value.  The  potential  range  using  this  threshold  value  includes  most 
of  the  state,  with  the  exception  of  the  Mojave  Desert,  the  very  north-central  and  northeast 
portion  of  the  state,  the  higher  elevations  of  the  Sierra  Nevada  Mountain  Range,  and  a 
southern  portion  of  the  San  Joaquin  Valley. 

The  Jenks  Natural  Breaks  method  produced  a  binary  map  with  suitable  habitat  above 
the  threshold  of  0.0530,  but  the  specificity  of  information  within  the  suitable  habitat  cat¬ 
egory  was  much  greater  than  in  either  of  the  previous  two  methods  (Tables  1,  2).  In  the 
map  classified  using  Jenks  Natural  Breaks  (Fig.  3),  the  two  highest  relative  habitat  cate¬ 
gory  rankings  (HCl  and  HC2)  are  centered  on  the  San  Francisco  Bay  area,  south  through 
the  Santa  Cruz  Mountains  to  habitats  on  the  Monterey  Bay  Peninsula,  east  through 
Sacramento  into  the  foothills  of  the  Sierra  Nevada  Mountains,  and  along  the  southern 
California  coasts.  Concentric  areas  of  decreasingly  suitable  habitat  (HC3  and  HC4)  sur¬ 
round  the  most  suitable  areas  and  continue  south  along  most  of  coastal  California  and  the 
coastal/inland  mountain  ranges  to  the  border  with  Mexico.  HC3  and  HC4  are  also  found 
along  river  corridors  exiting  the  Sierra  Nevada  Mountains.  Unsuitable  habitats  (HC5)  are 
in  and  east  of  the  Sierra  Nevada  Mountains,  the  northern  Cascade  Range,  a  major  portion 
of  the  San  Joaquin  Valley  south  of  Sacramento,  and  the  deserts  of  southern  California. 

Discussion 

Our  results  indicate  that  the  method  for  selecting  a  threshold  value  to  establish  a  break 
between  suitable  and  unsuitable  habitat  with  Maxent  is  an  especially  critical  issue  for  an  in¬ 
troduced  species  which  is  expanding  its  geographic  range.  We  have  shown  that  a  high  AUC 
value  by  itself  is  not  sufficient  to  support  the  predictive  accuracy  of  a  single  map  based 


RANGE  EXPANSION  OF  THE  EASTERN  GRAY  SQUIRREL 


63 


Uoi  I 


‘  A 


Future  Range  Estimate 

10th  Percentile  Training  Presence  Logistic  Threshold 

I  0  -  0.3051 
im  0.3051  - 1 

•  EGS  locations  2004  to  2018 

Sources:  Esri.  HERE,  Garmin,  Infermap,  increment  P  Corp.,  GEBCO,  USGS. 
FAO,  NPS,  NRCAN.  GeoBase,  IGN,  Kadaster  NL,  Ordnance  Survey.  Esri 
Japan,  METt.  Esri  China  (Hong  Kong),  swisstopo,  ©  OpenStreetMap 
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Fig.  1.  Estimate  of  suitable  habitat  predicted  by  the  model  using  the  10th  percentile  training  presence 
logistic  threshold  based  upon  eastern  gray  squirrel  locations  from  2004  to  2018. 

Upon  a  single  arbitrarily  selected  threshold  value.  AUC  is  not  a  perfect,  objective  measure 
of  the  predictive  power  of  the  model,  but  few  alternatives  are  available  for  presence  only 
data  (Merow  et  al.  2013).  The  commonly  used  10th  percentile  training  presence  logistic 
threshold  as  well  as  the  minimum  training  presence  logistic  threshold  are  clearly  arbitrary 
values  with  little  biological  basis  for  selection.  And,  as  shown  by  the  maps  presented  in  this 
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Fig.  2.  Estimate  of  suitable  habitat  predicted  by  the  model  using  the  minimum  training  presence  logistic 
threshold  based  upon  eastern  gray  squirrel  locations  from  2004  to  2018. 


paper,  these  thresholds  produce  vastly  dififereiit  projections  of  suitable  habitat.  The  use  of 
Jenks  Breaks  to  establish  a  threshold  value  is  based  upon  the  natural  clustering  of  values 
in  the  logistic  output  of  the  model  and  hence  should  be  less  arbitrary,  and  more  meaning¬ 
ful,  in  terms  of  the  biology  of  the  species  being  studied.  This  method  also  provides  graded 
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Table  I.  Threshold  breakpoints. 


10th  percentile  training 

Minimum  training  presence 

Jenks 

Class 

presence  logistic  threshold 

logistic  threshold 

Natural  Breaks 

Suitable 

0.3051  -  1 

.0299-1 

0.0530  -  1 

Unsuitable 

0-0.3051 

0-0.0299 

0  -  0.0530 

“likelihoods”  of  predicted  occupancy,  which  may  be  more  easily  tested  in  the  future  with 
new  occupancy  records. 

The  10th  percentile  training  presence  logistic  threshold  provided  the  most  conservative 
estimate  of  the  potential  range  by  assuming  that  some  of  the  presence  data  may  be  misiden- 
tified,  improperly  reported,  or  outside  of  the  area  in  which  the  EGS  can  persist,  and  then 
removing  those  points  (Young  et  al.  2011;  Uden  et  al.  2015).  The  map  is  inherently  skewed 
toward  the  realized  niche  of  the  EGS  in  California,  as  the  species  has  only  been  introduced 
to  a  few  locations  within  the  state.  It  is  closely  aligned  with  the  original  data  set  and  may 
be  biased  toward  human  accessible  areas.  The  threshold  may  be  too  conservative,  as  the 
EGS  is  still  spreading,  and  is  likely  to  tolerate  habitat  with  conditions  at  least  as  extreme 
as  those  in  the  current  range  in  California.  The  EGS  already  inhabits  areas  outside  of  the 
region  projected  as  suitable  in  this  map. 

The  map  created  with  the  minimum  training  presence  logistic  threshold  reflects  a  much 
broader  range  of  conditions  throughout  California  and  is  our  least  conservative  estimate  of 
the  future  range.  However,  use  of  this  threshold  value  most  likely  over-predicts  the  poten¬ 
tial  range  because  the  inclusion  of  a  single  erroneous  location  point  within  the  data  could 
greatly  affect  the  map.  The  map  using  this  lower  threshold  provides  an  accurate  estimate 
of  potential  range  only  if  it  is  absolutely  certain  that  all  location  data  have  been  correctly 
identified  as  EGSs,  and  all  individuals  and  populations  currently  exist  in  suitable  habitat. 
While  it  is  highly  likely  that  the  EGS  will  expand  its  range  to  include  some  of  the  areas 
mapped  by  this  method,  it  is  highly  unlikely  that  the  species  will  inhabit  all  of  the  area.  If 
the  EGS  is  introduced  to  or  expands  it  range  to  new  areas  and  survives,  the  distribution 
is  likely  to  more  closely  resemble  the  maximum  estimate  map  than  the  minimum  estimate 
using  the  10th  percentile  logistic  threshold. 

The  map  created  with  Jenks  Natural  Breaks  provides  the  most  useful,  detailed  classifica¬ 
tion  of  the  relative  habitat  suitability  rankings,  and  therefore  the  risk  of  invasion  to  various 
parts  of  California.  We  assume  that  the  lowest  ranking  habitat,  HC5,  is  unsuitable  based 
upon  the  lowest  grouping  of  values  (0.000  to  0.0530),  whereas  each  higher  rank  represents 
increasingly  suitable  habitat.  As  the  EGS  continues  to  spread,  the  habitat  suitability  rank¬ 
ing  may  increase  in  areas  with  habitat  similar  to  newly  invaded  areas  and  therefore  we 
cannot  say  that  this  map  is  a  permanent  ranking  of  habitat  suitability.  For  example,  areas 


Table  2.  Habitat  classes  created  with  Jenks  Natural  Breaks. 


Habitat  class 

Relative  habitat  suitability  ranking 

1 

0.4627  -0.7515 

2 

0.2917  -0.4627 

3 

0.1562-0.2917 

4 

0.0530-0.1562 

5 

0  -  0.0530 
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Fig.  3.  Estimate  of  suitable  habitat  predited  by  the  model  using  Jenks  Natural  Breaks,  with  five  habitat 
classes  (HC),  based  upon  eastern  gray  squirrel  locations  from  2004  to  2018.  HCl  represents  the  most  suit¬ 
able  habitat.  HC2-HC4  represent  decreasingly  suitable  habitat,  and  HC5  is  considered  unsuitable  habitat. 


now  ranked  as  HC3  could  in  the  future  be  ranked  as  HC2  if  many  populations  of  EGSs 
become  established  within  HC3.  We  believe  that  setting  the  threshold  for  suitable  habitat 
above  HC5  created  the  most  reasonable  estimate  of  the  potential  future  range  by  including 
a  broader  range  of  habitat  than  the  minimum  estimate,  but  excluding  the  most  extreme 
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desert,  mountain,  and  Central  Valley  habitats  of  the  maximum  estimate.  The  Jenks  Nat¬ 
ural  Breaks  map  is,  at  this  time,  our  most  reasonable  estimate  of  the  future  range  of  the 
EGS  in  California. 

We  selected  Maxent  to  create  the  ENM  because  of  its  high  performance  at  estimating 
local  occurrences.  Absence  data  were  not  available  or  informative,  so  a  presence-absence 
experiment  could  not  be  done  (Yackulic  et  al.  2012).  Additionally,  presence  only  data  may 
be  better  than  presence-absence  data  for  estimating  fundamental  niches  because  invasive 
species  have  not  yet  inhabited  their  full  potential  range  (Jimenez- Valverde  et  al.  2008).  For 
presence  only  data,  machine  learning  methods  such  as  Maxent  consistently  outperformed 
earlier  methods,  such  as  Bioclim  or  regression  models  in  predictive  success  (Elith  et  al. 
2006).  It  is  especially  good  for  small  and  medium  data  sets  because  the  generative  learning 
method  uses  an  algorithm  to  build  a  model,  as  opposed  to  a  discriminative  model  that 
estimates  the  value  of  one  categorical  variable  based  on  the  other  (Ng  and  Jordan  2001; 
Aguierre-Gutierrez  et  al.  2013). 

Location  data  were  limited  to  2004-2018  because  habitat  loss  has  been  extensive  in  Cal¬ 
ifornia  over  the  past  century,  to  ensure  the  use  of  current  species  distributions,  and  be¬ 
cause  satellite  collection  for  BIOCLIM  only  started  in  1972  (United  States  Geologic  Sur¬ 
vey  [USGS]  2015).  Overall,  2,825  of  the  3,725  EGS  records  in  California,  75.84%,  were 
obtained  from  wildlife  rehabilitation  records.  Records  indicated  the  location  that  the  per¬ 
son  submitting  the  squirrel  reported  finding  it.  Rehabilitation  facilities  are  very  likely  to 
accurately  identify  the  species,  but  may  cause  geographic  bias  toward  the  residential  areas 
surrounding  each  center.  Citizen  science  data  offers  the  benefit  of  having  broad  sources, 
but  may  be  skewed  toward  urbanized  areas,  roads,  or  other  easily  studied  sites  (Baldwin 
2009).  The  data  are  not  random,  but  rarefaction  has  been  used  to  equalize  the  sampling 
effort  among  areas,  which  allowed  us  to  make  a  reasonable  inference  of  the  species  distri¬ 
bution  (Yackulic  et  al.  2012).  The  assumption  that  detection  probability  is  constant  across 
sites  need  not  be  met  because  citizen  science  and  historic  records  simulate  repeat-visits  to 
each  site  (Yackulic  et  al.  2012). 

The  higher  elevations  of  the  Sierra  Nevada  were  classified  as  unsuitable  using  all  thresh¬ 
olds,  but  elevation  may  have  been  the  only  factor  causing  the  unsuitable  classification,  as 
the  highest  point  in  the  native  range  of  S.  carolinensis  is  at  6,643  feet,  at  Clingman’s  Dome 
in  Tennessee,  which  is  the  highest  point  in  the  state  and  the  Smoky  Mountains.  Elevation 
alone  may  not  actually  exclude  the  EGS.  The  species’  invasive  ability  worldwide,  especially 
into  areas  with  Mediterranean  climates,  such  as  South  Africa  and  Italy  (Gurnell  et  al.  2004; 
Bertolino  2008;  Benson  2013;  Bertolino  and  Lurz  2013)  suggests  that  it  could  acclimate  to 
conditions  in  much  of  California. 

Finally,  with  future  studies  the  potential  invasion  by  the  EGS  into  occupied  and  un¬ 
occupied  eastern  fox  squirrel  (EES)  or  WGS  habitat  is  an  excellent  test  case  for  several 
hypotheses  regarding  invasion  ecology.  These  hypotheses  include,  1)  “biotic  resistance”, 
which  suggests  that  high-biodiversity  ecosystems  are  more  resistant  to  invasion  than  low- 
diversity  ecosystems;  2)  “enemy  release”,  which  posits  that  the  absence  of  enemies  (i.e., 
competitors  and  predators)  increases  the  likelihood  of  invasion;  and  3)  “propagule  pres¬ 
sure/introduction  effort”,  which  proposes  that  the  introduced  population  size  and  the  fre¬ 
quency  of  introduction  can  contribute  to  successful  invasion  (Jeschke  2014).  The  first  hy¬ 
pothesis  could  be  tested  with  the  EES  by  comparing  the  rate  and  success  of  invasion  from 
neighboring  areas  of  low  and  high-biodiversity  systems  (e.g.,  oak  woodlands  in  the  Sierra 
Nevada  mountain  foothills).  The  second  hypothesis  could  be  tested  by  comparing  EGS 
invasion  and  persistence  success  with  the  presence  of  competitors  (e.g.,  EES  or  WGS) 
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and  predators,  such  as  hawks  and  owls.  Because  the  EGS  tolerates  urban  and  suburban 
areas,  predators  may  be  less  prevalent  in  areas  where  they  have  succeeded.  The  third  hy¬ 
pothesis  could  be  tested  using  a  combination  of  estimates  (e.g.,  from  historical  accounts) 
or  measurements  (e.g.,  for  newly-discovered  populations)  of  founder  population  size  and 
continued  introduction  (e.g.,  from  connected  populations)  and  reproductive  rate  to  predict 
expansion  and  persistence.  There  is  already  some  evidence  that  the  EGS  may  be  replac¬ 
ing  the  EES  on  the  western  side  of  the  San  Francisco  Bay  (Creley  and  Muchlinski  2017) 
so  assessing  potential  causal  factors  for  replacement  is  important.  Finally,  the  model  uses 
current  environmental  conditions,  and  could  be  adapted  to  include  estimates  for  changes 
in  climate,  which  may  reduce  the  amount  of  suitable  habit,  shift  it  to  higher  elevations,  or 
render  hotter  and  more  extreme  areas  of  currently  suitable  habitat  unsuitable  in  the  future. 
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The  San  Quintin  Kangaroo  Rat  is  Not  Extinct 


Scott  Tremor,^*  Sula  Vanderplank,’’^  and  Eric  Mellink^ 

^  San  Diego  Natural  History  Museum.  RO.  Box  121390.  San  Diego,  CA  92112 
^Departamento  de  Biologia  de  la  Conservacion,  Centro  de  Investigacion  Cientifica  y  de 
Educacion  Superior  de  Ensenada,  B.C.,  Carretera  Ensenada-Tijuana  #  3918,  22860 

Ensenada,  B.C.,  Mexico 


The  range  of  the  San  Quintin  kangaroo  rat  {Dipodomys  gravipes)  is  restricted,  so  far  as 
known  historically,  to  a  stretch  of  coastal  habitat  less  than  150  km  in  length  and  a  few 
kilometers  in  width  (Best  and  Lackey  1985)  at  the  southern  end  of  the  California  Floristic 
Province — a  global  biodiversity  hotspot  and  one  of  the  most  critically  endangered  ecosys¬ 
tems  on  earth  (Myers  et  al.  2000).  This  rodent  was  described  in  1925  by  Laurence  M.  Huey, 
who  reported  it  from  the  coastal  plains  from  San  Telmo  south  to  El  Socorro  and  on  the 
floodplain  of  the  Arroyo  El  Rosario  along  the  Pacific  coast  of  northern  Baja  California, 
Mexico.  The  area  between  El  Socorro  and  El  Rosario  is  largely  unsuitable. 

Huey  asserted  that  the  “mother  lode  of  this  species  is  found  near  Mesa  Agua  Chiquita, 
with  1,000  individuals  in  10  acres”  (field  notes  archived  at  the  San  Diego  Natural  History 
Museum).  He  described  the  habitat  in  this  area  as  hard  [clay]  soils  covered  with  grasses. 
Little  is  known  about  the  broader  habitat  requirements  for  D.  gravipes,  but  it  appears  sim¬ 
ilar  in  its  niche  requirements  to  Stephen’s  Kangaroo  rat  D.  stephensi  which  occupies  areas 
with  high  disturbance,  open  conditions  and  weedy  forbs  (Tremor  et  al.  2017).  Burrow  en¬ 
trances  were  closed  but  connected  by  visible  runways  up  to  75  m  long  (Fig.  1).  Previously, 
Nelson  (1922)  had  written  that  in  this  area  “the  vegetation  is  so  low  and  insignificant  that 
the  plain  has  the  appearance  of  an  open  prairie”.  Both  descriptions  could  refer  to  areas 
recovering  from  wheat  cultivation  years  earlier,  and  the  grasses  could  be  non-native  annual 
species  (e.g.,  Bronnis  spp.). 

Agriculture  in  the  San  Quintin  area  began  in  1891  when  British  farmers  converted  parts 
of  the  landscape  to  wheat  cultivation,  built  a  dam  for  irrigation,  and  installed  a  flour  mill 
(Taylor  1996).  The  settlement  was  abandoned  in  1917  (Phelts-Ramos  2004).  Subsequently, 
only  four  ranches  persisted  in  the  San  Quintin  Valley,  until  1947,  when  Title  3050  (which 
granted  agricultural  lands  to  families  from  other  regions  of  Mexico  as  a  cession  by  the  gov¬ 
ernment)  led  to  a  massive  expansion  of  agriculture  in  the  valley  (Ramirez- Velarde  2004). 
However,  full  expansion  of  agriculture  in  the  area  was  restricted  by  the  lack  of  roads  al¬ 
lowing  for  export  of  produce  from  the  area,  until  1973  when  the  road  connecting  it  with 
Ensenada  was  paved. 

In  1972,  before  the  building  of  the  transpeninsular  highway,  “the  broad  open  areas 
8.5  miles  N  of  San  Quintin  were  dotted  with  D.  gravipes  burrows”  (Best  1983).  Eight  years 
later  this  area  was  converted  to  cropland.  In  1980,  the  population  had  shrunk,  and  Best 
trapped  only  two  individuals  in  >1000  trap  nights.  Likewise,  the  area  9.6  km  east  of  El 
Rosario,  where  he  had  collected  35  specimens  in  1972,  was  in  1980  covered  by  the  paved 
transpeninsular  highway  (Best  1983).  The  Arroyo  del  Rosario  area  still  produced  7  speci¬ 
mens  in  1989  (Troy  L.  Best  in  litt.  to  EM,  20  July  1989). 


*  Corresponding  author:  stremor@sdnhm.org 
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Fig.  1.  Dipodowys  gravipes  habitat  at  Mesa  Agua  Chiquita.  Note  the  long  and  visible  runways.  Photo 
by  S.  Vanderplank. 


Dipodomys  gravipes  has  a  high  affinity  for  flat  terrain  and  is  intolerant  of  cultivation. 
Between  April  1989  and  September  1990,  Mellink  surveyed  all  habitats  where  the  species 
had  been  previously  found,  including  coastal  scrub,  fallow  fields,  and  river  wash.  Yet  this 
effort  totaling  more  than  800  trap-nights,  spread  over  12  cycles  of  trapping  events  within 
the  historical  range,  failed  to  yield  any  kangaroo  rats  other  than  D.  simukms.  As  a  result 
of  extensive  and  profound  habitat  alteration  by  agriculture,  D.  gravipes  has  been  listed 
as  endangered  by  the  Mexican  government  since  the  publication  of  its  first  list  of  species 
at  risk  (Instituto  Nacional  de  Ecologia  1994).  Since  then,  and  given  the  lack  of  further 
captures,  biologists  and  conservationists  have  feared  the  species  could  be  extinct  (Ceballos 
and  Navarro,  1991;  Mellink  1992',  1996-;  Mellink  and  Luevano  2005). 


*  Mellink,  E.  1992.  Status  de  los  heteromidos  y  cricetidos  endemicos  del  estado  de  Baja  California. 
Comunicaciones  Academicas.  CICESE.  Ensenada.  B.  C.  10pp. 

“Mellink,  E.  1996.  Problemas  de  conservacion  de  la  fauna  silvestre  en  el  estado  de  Baja  California. 
Conferencia  magistral.  XIV  Simposio  sobre  Fauna  Silvestre.  Facultad  de  Medicina  Veterinaria  y  Zootecnia. 
Universidad  Nacional  Autonoma  de  Mexico.  Mexico,  D.  F. 
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Fig.  2.  Dipodomys  gravipes  from  Mesa  Agiia  Chiqiiita.  Note  the  rear  foot  measurement  of  46  mm  and 
thin  white  lateral  tail  stripe.  Photo  by  S.  Vanderplank. 


On  4  July  2017,  Tremor  and  Vanderplank  during  routine  trapping  and  general  inven¬ 
tory,  placed  traps  in  the  vicinity  of  San  Telmo  and  the  vicinity  of  San  Quintin.  Sherman 
traps  (30.5  cm)  were  set  from  sundown  (approximately  20:00  PM)  to  dawn  (approximately 
05:30  AM)  and  baited  with  rolled  oats.  They  set  13  traps  along  runways  and  burrow  en¬ 
trances  on  a  disturbed  embankment  adjacent  to  a  fallow  agricultural  field  5.6  km  east  of 
San  Quintin,  on  Mesa  Agua  Chiquita  (near  the  “mother  lode”  described  by  Huey).  Runs 
were  visible  in  the  non-native  vegetation,  which  included  Meseinhiycmthemum  crystalliniim, 
Brassica  tournefortii,  Salsola  tragus,  and  Hirschfeldia  incana,  with  individuals  of  the  native 
species  Ambrosia  chenopodiifolia  and  Malosma  laurina  nearby.  Three  female  and  one  male 
D.  gravipes,  all  adults,  were  trapped.  Other  species  trapped  were  three  individuals  of  the 
pocket  mouse  Chaetodipus  fallax  and  one  woodrat  {Neotoma  hryanti).  The  pocket  mice  ap¬ 
peared  to  be  sharing  the  burrows  of  D.  gravipes.  Traps  were  set  at  three  additional  nearby 
locations  and  at  San  Telmo,  but  no  D.  gravipes  were  caught  there. 

Identification  of  D.  gravipes  was  based,  among  other  characters,  on  measurements  of 
females,  greater  than  the  range  of  values  of  D.  sinndans  (Table  1).  Five  toes  on  its  hind  feet 
distinguish  Dipodomys  gravipes  from  sympatric  D.  merriami,  which  have  four  toes.  Also, 
diagnostic  was  the  width  and  pattern  of  the  tail  which  in  D.  gravipes  is  thicker  than  in 
D.  simukms,  with  the  white  lateral  tail  stripe  narrow  and  indistinct  where  it  merges  with 
black  dorsal  and  ventral  tail  stripes  (Best  and  Lackey  1985).  The  feet  were  also  thick,  as  de¬ 
scribed  by  Huey  (1925).  The  adults  captured  were  unusually  feisty  for  the  genus  Dipodomys 
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Table  1.  Body  measurements  (mm)  of  D.  gravipes  females  captured  (2017)  and  females  only  specimen 
data  from  the  San  Diego  Natural  History  Museum  (SDNHM)  from  the  same  species  and  its  closest  con¬ 
generic  species.  D.  siiiiukms.  All  measurements  in  millimeters.  Ear  length  excluded  as  historic  methods  of 
taking  these  measurements  ditfers  from  current  methods. 


Dipodomys  gravipes  D.  simukms 


SDNHM  SDNHM 


Three  captures  2017 

specimens  from 

San  Quintin 
region 
(n  =  35) 

Best  1983 
Mean 

specimens  from 
San  Quintin 
region 
(n  =  51) 

1 

2 

3 

Mean 

Total  length 

295 

290 

300 

295 

280-320 

300 

247-292 

Tail 

175 

185 

188 

183 

161-187 

173 

145-175 

Hind  foot 

46 

42 

40 

43 

41-45 

44.1 

38-42 

and  difficult  to  handle;  the  adult  male  escaped  before  being  measured,  though  it  appeared 
significantly  larger  than  the  females. 

Data  from  the  live  captures  were  tabulated  and  later  compared  to  specimens  housed  at 
the  San  Diego  Natural  History  Museum  (SDNHM),  which  houses  extensive  collections 
of  this  species.  Museum  collections  of  D.  sinndans  from  the  local  region  were  compared 
to  D.  gravipes  in  order  to  provide  a  more  accurate  comparison,  since  D.  simukms  has  a 
wide  range  and  can  be  highly  variable.  Ear  measurements  were  not  included  because  they 
have  historically  been  recorded  using  differing  techniques  and  therefore  cannot  be  used  for 
accurate  comparison. 

The  presence  of  this  species  adjacent  to  fallow  agricultural  land  suggests  that  the  recent 
drought  and  subsequent  fallowing  may  have  increased  available  habitat  and  favored  the 
recovery  of  D.  gravipes.  If  this  is  true,  concern  should  be  elevated  for  this  species  in  the 
event  of  increased  intensity  of  agriculture  in  future,  if  and  when  the  drought  abates.  It  is 
not  clear  whether  this  discovery  represents  a  rebound  of  this  species  or  a  last  remnant 
population.  Surveys  have  begun  throughout  the  broader  region  in  an  attempt  to  guide 
future  efforts  to  recover  the  species. 

This  discovery  is  a  note  of  hope  for  the  species  and  of  great  importance  to  regional 
conservation.  Nevertheless,  the  primary  threats  of  habitat  destruction  and  agricultural  ex- 
pansion  (Alvarez-Castaheda  et  al  2008)- ,  which  cause  fragmentation  and  reduced  genetic 
vigor,  in  addition  to  direct  habitat  loss,  leave  this  species  highly  threatened.  The  human 
population  in  the  range  of  D.  gravipes  also  continues  to  increase,  presumably  increasing 
pressures  such  as  artificial  lighting,  application  of  rodenticides,  predation  by  domestic  pets 
and  overgrazing  by  livestock.  Additional  data  on  range  and  abundance  are  needed,  and 
genetic  studies  in  areas  previously  occupied  are  recommended  in  the  absence  of  voucher 
specimens.  Our  finding  is  undeniably  good  news,  but  no  conservation  strategy  can  be  pro¬ 
posed  until  the  current  status  and  distribution  of  D.  gravipes  is  better  understood. 


Alvarez-Castaiieda,  S.T.,  Castro-Arellano,  I.  &  Lacher,  T.  2008.  Dipodomys  gravipes.  The 
lUCN  Red  List  of  Threatened  Species  2008:  e.T6676A  12794061.  Downloaded  on  02  October 
20 17.http://dx.doi.org/10.2305/IUCN.UK.2008.RLTS.T6676A  1279406  Ten. 
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The  lowlands  of  cismontane  southern  California  have  proven  to  be  hospitable  not  only  to 
humans  but  to  many  exotic  plant  and  animal  species  (Cox  1999).  In  addition  to  14  nonna¬ 
tive  reptile  species  established  in  the  area'  is  a  localized  population  of  confusing  whiptails 
in  Orange  County  that  was  first  reported  by  Winkleman  and  Backlin  (2016).  That  report 
was  based  on  observations  in  south  Irvine  in  May  2014  and  April-June  2015  and  at  least 
one  similar  whiptail  seen  in  adjacent  Lake  Forest  in  July  2015.  Four  specimens  collected  in 
Irvine  at  that  time  were  identified  as  belonging  to  the  Aspidoscelis  flagellicauda/ sonorae 
complex.  The  Gila  Spotted  Whiptail  {A.  flagellicauda)  and  Sonoran  Spotted  Whiptail 
{A.  sonorae)  are  morphologically  similar  all-female  species  native  to  Arizona,  New  Mexico, 
Sonora,  and  Chihuahua  and  previously  not  known  to  occur  away  from  their  native  ranges. 

Gary  Nafis“  provided  an  update  on  the  status  of  these  lizards,  noting  that  they  are  not 
extirpated  in  Irvine,  as  suspected  by  Winkleman  and  Backlin  (2016),  and  “as  of  7/17,  they 
have  been  found  only  in  Orange  County  in  Irvine,  Lake  Forest,  and  Aliso  Viejo,  but  they 
appear  to  be  spreading  quickly.”  The  reference  to  Aliso  Viejo  was  presumably  based  on  the 
observations  detailed  here.  Information  displayed  on  iNaturalist'^  under  the  name  of  Sono¬ 
ran  Spotted  Whiptail  includes  the  Irvine  observations  and  others  beginning  in  2015,  but 
the  precise  locality  data  is  “obscured”  according  to  the  wishes  of  reporting  individuals  or 
institutions.  All  of  the  iNaturalist  locations,  obscured  or  otherwise,  are  within  the  general 
geographic  boundaries  described  above  (G.B.  Pauly  pers.  comm.). 

Our  observations  come  primarily  from  a  church  campus  in  suburban  central  Laguna 
Woods  (33.609882  N,  -1 17.733124  W;  adjacent  to  the  vast  retirement  community  formerly 
known  as  Leisure  World),  approximately  4.8  km  south-southwest  of  the  locations  reported 
by  Winkleman  and  Backlin  (2016).  The  manicured  landscape  with  scattered  ornamental 
shrubs  provides  suitable  habitat  for  the  whiptails,  not  unlike  the  situation  described  by 
Winkleman  and  Backlin  (2016).  Sandy  substrates  are  especially  favored  by  these  lizards. 
Like  Winkleman  and  Backlin,  we  initially  struggled  with  their  identification,  generally  try¬ 
ing  to  force  the  label  of  Orange-throated  Whiptail  {A.  liyperythra)  upon  them. 

Our  first  observation  was  of  a  single  individual  on  13  June  2010.  Infrequent  observations 
continued  through  2016,  but  in  2017  we  increased  our  effbrt  to  document  them.  Spotted 


*  Corresponding  author:  richard.erickson@lsa.net 
http://www.californiaherps.com  (accessed  May  2018). 
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'  Aspidoscelis  fiagellicauda  x  Aspidoscelis  sonorae  complex,  http://www.californiaherps.com  (accessed 
May  2018). 
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Fig.  1 .  Spotted  Wliiptail  {Aspidoscelis  fiageUicauda/souorae  complex)  in  residential  Aliso  Viejo,  Orange 
County,  1  June  2018.  Photo  by  Richard  A.  Erickson. 


whiptails  were  seen  from  16  April  2017  to  22  July  2017  but  were  not  seen  on  30  July  2017 
or  thereafter.  We  collected  individuals  on  1  and  4  July  2017  (LACM  189582  &  189583; 
collected  under  California  Department  of  Fish  and  Wildlife  Scientific  Collecting  Permit 
#000777)  and  saw  at  least  three  more  after  that.  Our  first  observation  of  2018  was  on 
22  April  2018. 

On  three  occasions  in  May  and  June  2017,  and  three  more  in  May  and  June  2018  (Fig. 
1),  we  observed  spotted  whiptails  amongst  ornamental  landscaping  at  two  locations  in 
nearby  portions  of  Aliso  Viejo,  along  Canyon  Wren  Lane  and  Calle  Corta  to  the  inter¬ 
section  of  Chickadee  Lane.  These  observations  extended  the  known  range  another  1 .2  km 
south-southwest  of  the  church  site. 

In  summary,  all  of  our  observations  of  spotted  whiptails  were  in  suburban  landscaping 
>0.6  km  removed  from  native  coastal  sage  scrub  and  up  to  6  km  south-southwest  of  the 
original  Irvine  location.  Our  observations  of  activity  extend  from  16  April  2017  to  22  July 
2017  with  the  exception  of  one  fall  sighting  on  4  September  2016  (iNaturalist  records  ex¬ 
tend  to  31  October  2017).  In  light  of  the  apparent  displacement  of  Western  Fence  Lizards 
(Sceloporus  occidentalis)  by  nonnative  Italian  Wall  Lizards  {Podavcis  sicidiis)  elsewhere  in 
in  the  Los  Angeles  Basin  (G.B.  Pauly  pers.  comm.)  and  the  related  warning  of  Deichsel 
et  al.  (2010),  we  note  that  as  of  May  2018  fence  lizards  were  still  present  at  all  locations 
where  we  have  seen  whiptails.  The  spotted  whiptail  population  in  Orange  County  is  known 
to  have  been  present  since  at  least  2010,  with  the  first  observations  in  the  city  of  Laguna 
Woods.  It  has  since  been  found  in  the  neighboring  cities  of  Irvine,  Lake  Forest,  and  Aliso 
Viejo.  This  area  has  received  industrial-scale  residential  landscaping  since  the  1960s.  The 
widespread  movement  of  plants,  soil,  and  materials  involved  in  such  endeavors  has  been 
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implicated  in  the  establishment  in  southern  California  of  species  such  as  the  Brown  Anole 
{Aiiolis  sagrei;  Mahrdt  et  al,  2014),  Common  Coqui  [Eleutherodactylus  coqui^;  SDSNH 
76135  and  76138),  and  others  (G.B.  Pauly  pers.  comm.).  While  we  hope  that  these  species 
will  remain  restricted  to  the  highly  altered  urban  environment,  the  spotted  whiptail  has  al¬ 
ready  been  found  at  the  edge  of  protected  open  space  in  Irvine,  where  it  is  feared  that  it  may 
threaten  or  outcompete  native  species  in  natural  habitats  (R.S.  Winkleman  pers.  comm.). 
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